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(II.)  Disclaimer 

The  findings  in  this  report  are  solely  those  of  the  author 
and  are  not  to  be  interpreted  as  the  official  position  of  the 
Air  Force  of  Scientific  Research  or  the  U.  S.  Government. 
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This  document  constitutes  the  report  on  Grant  AFOSR-81 -01 70 
for  the  period  from  March  15,  1981  to  March  14,  1982.  The  format 
under  which  the  report  is  organized  is  described  as  follows. 

There  were  three  major  and  separate  research  tasks  carried  out 
during  the  past  year  under  the  sponsorship  of  Grant  AFOSR-81 -01 70 . 
These  separate  research  task  areas  were:  optical  sub-sampling 
for  mul ti spectral  image  data  compression;  opti cal / d  i  gi ta 1  hybrid 
system  architectures  for  interframe  (temporal)  image  data  compres¬ 
sion;  the  usage  of  radiometric  and  spatial  trans forma ti ons  for 
adaptive  image  data  compression.  Since  each  of  these  research 
task  areas  is  independent  of  the  others,  we  have  chosen  a  report 
format  which  emphasizes  the  independence  and  makes  it  possible 
to  concentrate  upon  reading  about  the  activity  in  a  particular 
topic  area  with  no  involvement  or  conflict  with  topical  material 
from  another  area.  Thus,  the  major  research  results  are  presented 
in  Appendices  I,  II,  and  III,  three  appendices  that  can  be 
separated  from  the  rest  of  the  report,  if  so  desired.  As  such, 
each  appendix  could  be  considered  a  "mini -report" .  Each  ap¬ 
pendix  has  its  own  introductory  section  which  summarizes  the 
research  objectives  of  the  research  reported  in  that  appendix. 
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Optical  Sub-Sampling  for 
Mul ti spectral  Image  Data  Compression 
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Optical  Sub-Sampling  for  Mul ti spectral 
Image  Data  Compression 


Research  Objectives 

Conventional  methods  in  mul ti -spectral  image  data  compression 
require  de-correl ation  trans format! ons ,  e.g.,  the  Karhunen-Loeve, 
which  are  very  costly  in  computing.  The  objective  of  this  task  is 
to  demonstrate  that  substantial  data  compression  (e.g.,  4:1)  can 
be  achieved  with  little  sacrifice  in  quality.  The  method  is  to 
under  sample  each  image  array  of  a  mul ti spectral  system,  but  with 
the  sample  array  of  each  spectral  band  staggered  or  offset  so  that 
sample  pixels  missing  in  a  given  band  are  sampled  in  other  bands. 
The  use  of  the  combined  spatial  and  spectral  correlation  within  and 
between  bands  is  used  to  estimate  missing  pixels  and  reconstruct 
imagery  of  higher  quality  than  that  obtained  solely  from  under¬ 
sampled  data. 

Appendix  Contents: 

This  appendix  consists  of  a  summary  of  the  research  in  this 


area  as  directly  evaluated  to  date,  plus  a  paper  on  general 
problems  in  image  sampling  which  is  cf  relevance  to  the  sub¬ 
sampling  which  is  employed  in  this  method. 
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This  section  of  the  report  will  review  progress  on  the  project 
to  achieve  data  compression  of  mul ti spectral  imagery  by  staggered 
array  sensor  undersampling  and  subsequent  reconstruction  of  the 
original  scene  utilizing  interband  redundancy  of  edge  i nf orma ti on . 
The  technique  assumes  that  features  possessing  high  spatial  fre¬ 
quencies  are  similar  in  all  spectral  bands  of  a  scene.  This  will 
be  true  for  shadow  borders;  however,  spectrally  dependent  edges 
such  as  the  junction  of  materially  different  terrain  regions  may 
be  less  spectrally  correlated. 

Prototype  solid-state  cameras  have  been  constructed  which  use 
similar  staggered  color  (RGB)  arrays  [1],  [2] ,  [3],  The  device 
in  [1]  comprises  a  set  of  adjoining  linear  CCD's,  each  array 
being  sensitive  to  only  one  spectral  band.  The  arrays  are  de¬ 
ployed  in  alternating  color  sensitivity,  with  elements  of  each 
strip  offset  from  the  elements  of  neighboring  strips.  Dye-deposi¬ 
tion  color  filter  arrays  (CFA's)  have  been  bonded  to  sensor  array 
chips  [2],  13].  The  CFA's  are  rectilinear  masks  with  repeating 

pattern  R  G;  the  mask  is  thus  staggered  by  spectral  element  but 
G  B 

without  sensor  gaps  as  was  the  case  in  [1].  Configurations  up  to 
484  x  384  elements  with  34  pm  x  20  pm  element  size  have  been 
achieved  [3],  However,  little  attempt  has  been  made  to  estimate 
imagery  in  unsampled  bands  from  data  in  a  sampled  band:  in  [2], 
edge  information  in  the  high  resolution  G  band  was  added  directly 
to  the  low  resolution  P.  and  3  bands. 

The  sampling  scheme  considered  in  the  present  research  con¬ 
sists  of  four  bands  of  a  mul ti s pectra 1  sensor  which  are  equally 
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undersampl ed ,  although  at  full  resolution,  at  staggered  intervals 
in  a  scene.  The  result  is  a  single  composite  image,  a  mosaic  of 
the  band  samples,  having  identical  spatial  dimensions  to  the  original 
scene.  The  four  spectral  estimates  of  the  scene  are  then  to  be  re¬ 
constructed  in  full  detail  from  this  single  image.  The  spatial 

arrangement  of  pixels  in  a  four-band  image,  bands  I  -  IV,  is 
» 

I  II  I  II 

III  IV  III  IV 

I  II  I  II 

III  IV  III  IV 


A  reconstruction  window  five  samples  square  is  passed  across  the 
mosaic  image;  at  each  window  location,  the  center  sample,  plus 
estimates  of  the  other  three  bands  formed  from  the  window  content, 
are  written  to  the  correspondi ng  spatial  coordinates  in  four 
output  images. 

The  contents  of  the  window  may  be  described  by  the  array: 

X  H  S  H  X 

V  0  V  D  V 

S  H  S  H  S 

V  D  V  D  V 

X  H  S  H  X 

where 

S  is  the  band  of  the  central  sample, 

H  is  the  band  of  the  horizontally  nearest-neighbor  sample, 


V  is  the  band  of  the  vertically  nearest-neighbor  sample, 

D  is  the  band  of  the  diagonally  nearest-neighbor  sample,  and 

X  positions  are  ignored. 

In  each  window  position  the  value  of  S  is  given,  and  the 
values  of  H,  V,  and  D  must  be  estimated  from  the  contents  of 
the  window.  If  E  is  the  general  band  to  be  estimated,  ( i . e .  ,  H,  V, 
or  D)  then  the  values  S  and  E  may  be  decomposed  into  local  means  and 
differences  as: 

S  =  S  +  AS 

E  =  E  +  AE.  (1) 

A 

Let  us  form  the  estimate  E  =  F  +  aAS,  i.e.,  we  wish  to  estimate 

*  Q 

AE  from  AS.  A  MMSE  criterion  such  that  <(e  -  E)  >  a  minimum  will 
be  assumed.  The  band  images  of  E  and  S  will  be  treated  as  weakly 
stationary  random  variables  over  the  extent  of  the  estimation 
window,  so  that  the  ensemble  averages  are  replaced  by  spatial 
averages.  The  coupling  parameter  a  resulting  from  the  MMSE 
optimization  is 

a  =  PES(oE/as),  (2) 

where  a£  and  cr^  are  the  standard  deviations  of  the  E  and  S 
bands,  respectively,  and  is  the  correlation  coefficient 
between  E  and  S,  defined  by 


'ES 


<AEAS> 
aE  aS 


”  1  — p  E  S—  1  • 


(3) 


Note  that  negative  values  of  o  account  for  contrast  reversals 
between  bands,  and  that  for  c  =  0  (uncorrelated  E  and  S  bands) 
the  band  estimate  is  the  local  mean.  If  a  is  correctly  chosen 
according  to  (2),  then  the  residual  mean  square  error  depends 
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on  PES: 

=  oE2(1-pES2)  .  (4) 

When  pE^  =  -  1,  perfect  reconstruction  is  possible;  when  E  and  S 

o 

are  uncorrelated  (and  when  a  =  0),  then  t  equals  the  variance 
of  the  estimated  band.  The  typical  situation  will  fall  between 
these  extremes. 

The  estimation  of  the  correlation  coefficient  pES  at 
each  window  location  is  non-trivial.  If  the  E  and  S  band  values 
of  each  window  element  were  a  priori  known  (in  which  case  no 
recons tructi on  would  be  needed),  the  expectation  <aEaS>  could  be 
formed  from  the  average  product  of  the  E  and  S  differences  at 
each  element  due  to  the  assumption  of  local  stati onari ty .  In 
the  absence  of  such  full  information,  an  initial  approximation 
of  pES  may  be  obtained  by  averaging  the  product  of  interpolated 
values  at  each  window  element.  The  accuracy  of  the  pe^  approxima 
tion  will  then  depend  on  the  within-band  spatial  correlations  of 
both  E  and  S. 

Simulation 

The  above  algorithm  has  been  used  in  the  reconstruction 

of  a  demonstration  image.  The  upper  left  photograph  in  figure  1 

shows  the  original  green  band  of  a  tricolor  (RGB)  256x256  8-bit 

digitized  image.  In  place-  of  four  distinct  bands,  the  mosaic  will 

be  formed  from  a  sampling  array  R  G  ,  with  the  R  red  band  in 

B  R 

both  positions  I  and  IV.  The  three  associated  c  maps  (G-Rj, 
G-Rjy,  G-B)  are  included  in  figure  1  clockwise  from  the  top.  A 
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closeup  of  the  original  is  given  in  figure  2,  and  figure  3 
displays  the  mosaic  sensor  image. 

It  is  often  true  that  some  pairs  of  bands  are  more  spectrally 
correlated  than  others.  This  knowledge  can  be  used  in  forming 
the  final  band  estimate.  For  example,  in  the  demonstration  image 
bands  G-B  and  bands  Rj-Rjy  are  such  pairs.  A  post-reconstruction 
filter  can  be  devised  that  replaces  an  estimate  by  its  within- 
band  nearest  neighbors  average  if  the  local  is  below  a 

threshold  and  if  the  local  S  and  E  bands  are  not  hi gh-correlation 
pairs.  The  sample  array  symmetry  then  guarantees  that  the  nearest- 
neighbor  within-band  estimates  were  formed  from  high-correlation 
pairs  and  tend  to  be  more  reliable  estimates.  Thus  the  best 
surrounding  reconstruction  estimates  will  supplant  a  poor  estimate. 

Figure  4  shows  the  result  of  reconstruction  of  the  mosaic- 
sampled  image  of  figure  2,  after  the  described  pos t- f i 1 teri ng 
with  a  | p g ^ i  threshold  of  0.5.  Approximations  of  pS£  were 
generated  from  b i 1 i near- i n terpol a  ted  versions  of  the  mosaic.  For 
comparison,  a  direct  bilinear-interpolation  of  the  mosaic  sampled 
original  displays  considerable  blurring  (figure  5).  Cubic 
spline  and  sine  interpolation  results  are  similarly  degraded. 

The  root  normalized  mean  square  error  (RNMSE)  between  the 
original  and  estimate  of  each  band  is  given  in  the  following 


table: 
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Demonstration  Image  -  RNMSE 

R-I 

G-II 

B-III 

R-IV 

Reconstruction , 

"true  p"  maps 

.0291  7 

.06127 

.06590 

.02869 

Post- f i 1 tered  version  of 

above 

.02822 

.05960 

.06464 

.02773 

Reconstruct!' on  , 

Interpolated  p  maps 

.02993 

.06235 

.06941 

.02931 

Post-filtered  version 

of  above 

.02875 

.06161 

.06842 

.02832 

Bilinear  interpolation 

.03247 

.07286 

.07605 

.03206 

The  correlati on-reconstructions  are  seen  to  be  more  accurate 

than 

the  bilinear  interpolation. 

The  replacement  post-filtering 

with  a  |  p  g  F  !  threshold  of  0. 

5,  makes  < 

j  small  improvement  in 

the 

estimation  of  all  bands. 


$  umma  ry 

Test  results  have  indicated  that  interband  correlation  is 
useable  in  a  simple  mul ti spectral  data  compression  scheme  whose 
only  coding  is  the  spatial  configuration  of  the  sensor  array. 

The  4:1  maximum  compression  ratio  for  the  four-band  case  does 
not  take  into  account  any  further  compression  available  by  the 
usual  single-band  methods,  e.g.,  DPCM.  Future  research  will 
consider  the  effect  of  the  sampling  scheme  on  single-band  methods. 
The  performance  of  various  recons tructi on-wi ndow  configurations, 
the  treatment  of  ncnstati onari ty  within  the  window,  and  simula¬ 
tions  on  a  more  varied  selection  of  imagery  will  also  be  examined. 


Figure  1.  Green  Sand  Original  with  Correlation  Coefficient  Maps 


Figure  L  . 


Detail  of  Original,  Single-Band 
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Abstract 

This  paper  is  a  one-dimensional  analysis  of  the  degradation  caused  by 
image  sampling  and  interpolative  reconstruction.  The  analysis  includes 
the  sample-scene  phase  as  an  explicit,  random  parameter  and  provides  a 
complete  characterization  of  this  image  degradation  as  the  sum  of  two  terms; 
one  term  accounts  for  the  mean  effect  of  undersampling  (aliasing)  and  non-ideal 
reconstruction  averaged  over  all  sample-scene  phases;  the  other  term  accounts 
for  variations  about  this  mean.  The  results  of  this  paper  have  application 
to  the  design  and  performance  analysis  of  image  scanning,  sampling  and 
reconstruction  systems. 
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Introduction 

In  1954  Mertc  and  Gray^  published  their  pioneering  paper  on  the  performance 
of  a  particular  line-scan  imaging  system,  television.  Over  the  next  20  years 
this  work  was  extended  considerably  by  Schade"  and  others  ’  .  More  recently 
a  diversity  of  line-scan  imaging  systems  have  evolved  from  the  application 
of  digital  and  electro-optical  technology  in  such  areas  as  space  exploration, 
remote  sensing,  reprographics,  and  medicine.  These  more  recent  scanning 
systems  typically  employ  electronic  sampling  and  digital  (or  optical)  image 
reconstruction  and  their  performance  has  been  analyzed  in  numerous  papers 
and  texts,  for  example  references  5  through  12.  Although  various  effects 
have  been  considered  in  these  performance  analyses,  one  which  is  commonly 
ignored  is  that  associated  with  the  phase  uncertainty  of  the  scene  relative 
to  the  sampling  grid. 

It  is  well  known  that  for  practical  systems,  image  sampling  and  reconstruc- 

1 

tion  inevitably  produce  seme  degradation  in  the  reconstructed  image  ;  the 

sampling  contribution  to  this  degradation  is  the  Moire-like  phenomenon  known 

as  aliasing  and  the  reconstruction  contribution  is  a  blurring  caused  by 

high  frequency  supression.  Aliasing  has  been  investigated  for  line-scan 

13  14 

systems  both  experimentally  and  theoretically  ,  but  the  effect  of  sample- 
scene  phase  was  not  addressed.  Similarly  this  phase  effect  is  typically 
ignored  in  studies  of  various  reconstruction  filters,  e.g.,  reference  15.  The 
amount  and  appearance  of  the  image  degradation  associated  with  sampling  and 
reconstruction  will  change  as  the  sampling  grid  is  shifted  across  the  scene, 
i.e.,  as  the  sample-scene  phase  is  varied.  How  does  the  magnitude  of  this 
degradation  depend  upon  the  sample-scene  phase?  To  date  a  complete  analysis 
of  this  question  has  not  appeared. 
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In  this  paper  the  one-dimensional  process  of  image  scanning,  sampling, 
and  interpc-lative  reconstruction  is  analyzed  in  both  the  spatial  and  frequency 
(Fourier)  domain  with  the  sample-scene  phase  appearing  as  an  explicit, 
random  parameter.  To  produce  mathematically  tractable  results  in  the  spatial 
domain  it  is  necessary  to  restrict  attention  to  simplified  models  of  the 
scene,  scanner  PSF,  and  reconstruction  function.  Nevertheless ,  considerable 
insight  into  the  general  problem  is  gained  with  this  spatial  domain  analysis. 
The  results  obtained  in  the  frequency  domain  are  completely  general 
with  no  restriction  on  the  scene  frequency  spectrum,  scanner  MTF,  or  recon¬ 
struction  filter.  This  approach  character! :es  completely  the  relation  of  the 
sample-scene  phase  to  scanning,  sampling,  and  reconstruction  and  provides 
several  measures  of  the  resulting  image  degradation  which  can  be  easily 
evaluated  by  numerical  quadrature.  It  is  demonstrated  that  although  image 
degradation  due  to  sampling  and  reconstruction  is  present  for  all  practical 
systems,  its  average  magnitude  can  be  controlled  by  a  proper  choice  of  scanner 
MTF,  sampling  rate,  and  reconstruction  filter. 

Although  the  results  of  this  paper  have  general  application,  the  primary 
motivation  for  the  research  was  t«>  incorporate  the  effect  of  sample-scene 
phase  into  an  analysis  of  the  image  degradation  due  to  sampling  and  recon¬ 
struction  of  remotely-sensed  multispectral  (e.g.,  Landsat)  images.  For  such 
images  there  is  known  to  be  (see,  for  example,  reference  16)  a  sample-scene 
phase-dependent  blurring  and  apparent  location  shift  of  high  contrast 
features  due  to  sampling  and  reconstruction.  When  classifying  multispectral 
image  data,  one  is  concerned  with  any  blurring  and  shift  in  location  of  edges, 
lines,  points,  etc.,  since  such  high  contrast  features  are  commonly  used  as 

i  _ _ 
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control  points  in  the  resampling  process  of  registering  one  image  to  another, 
and  even  small  shifts  in  control  point  locations  can  result  in  significant 
misregistration  and  subsequent  misclassifications * ' . 

Formulation 

Figure  (la)  illustrates  a  one-dimensional  scene,  denoted  f(x-u),  which 
contains  a  high  contrast  feature,  in  this  case  an  edge,  at  the  point  u.  A 
typical  scanner  ?SF,  denoted  h(x),  is  illustrated  in  Fig. (lb)  and  the  image 

g(x-u)  =  h(x)*f(x-u)  (1) 


formed  by  convolving  the  scene  with  the  PSF  is  illustrated  in  Fig.  (ic). 
The  process  of  image  sampling  can  be  represented  symbolically  as  the 
product 


g(x-u)  comb(x)  (2a) 

18  19 

where  sampling  is  accomplished  by  the  familiar  comb  (or  Shah  )  function 


comb  (x)  *  £  <5  (x-n)  . 

n=-=° 


(2b) 


* 

The  spatial  coordinate  x  is  normalised  in  units  of  sample  interval  so  that 
image  sampling  occurs  at  the  integer  locations  3S  indicated  in  Fig.  (Id).  In 
this  coordinate  system  u  is  the  sample-scene  phase  parameter;  values  of  u 
between  zero  and  one  indicate  the  position  (or  phase)  of  the  scene  relative 
to  the  sampling  grid. 

12 

Image  reconstruction  is  commonly  modeled  as 


gr(x;u)  =  [g(x-u)  comb(x)]  *r(x) 


(3a) 
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where  the  interpolation  function  r(x)  is  the  impulse  response  of  the  (linear 
shift  invariant)  reconstruction  process.  A  typical  interpolation  function 
and  the  resultant  reconstructed  image  gr(x;u)  are  illustrated  in  Figs,  (le) 
and  (lf)>  respectively.  Equation  [3a]  can  be  written  equivalently  as 


SrC*;u)  =  ]T  g(n-u)r(x-n) 

n*-»  (oo) 

which  demonstrates,  consistent  with  the  notation  g„(x;u),  that  the  reconstructed 

tmage  is  not  simply  a  function  of  the  difference  x-u.  Equation  C3b3  also 

illustrates  (along  with  Fig.  (If))  that,  for  a  fixed  x,  the  extent  to 

which  neighboring  image  samples,  g(n-u),  contribute  to  the  reconstructed 

image  is  determined  by  the  spread  cf  the  interpolation  function.  For  digital 

image  reconstruction  (sometimes  called  resampling”  ),  computational 

considerations  dictate  that  this  spread  should  be  small;  typically  r(x)  is 

identically  :ero  for  :x|  larger  than  2  or  5. 

Scanning  produces  an  image  (Fig. (lc))  which  is  a  blurred  (i.e.  smoothed) 
copy  of  the  original  scene  (Fig.  (la));  similarly,  subsequent  sampling  and 
reconstruction  produces  additional  blurring  in  the  reconstructed  image 
(Fig-  (!*))■  to  study  the  blur  (i.e.,  the  square  of  the  radiometric  error) 

present  in  a  reconstructed  image  it  is  convenient  to  define  two  quantities 


image  blur  »  -c~ 

and 


>30 


(4) 


sampling  and  reconstruction  blur 


These  two  quantities  are  associated  with  (but  not  equal  to)  the  shaded  areas 
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indicated  in  Figs.  (2a)  and  (2b);  and  s|R  represent  the  energy  of 
the  differences  f (x-u) -g(x-u) ,  and  g(x-u)-gr(x;u) ,  respectively.  How  do 

-cn  -depend  upon  the  sample-scene  phase  parameter  u?  For  s?  the  answer  is 


'I 


obvious,  for  s|R,  it  is  not. 


Image  £lur  (Eq.  [4])  is  an  inevitable  but  well  understood  effect  caused 
solely  by  scanning  the  scene  with  a  non-ideal  (i.e.,  non-impulsive)  FSF.  Image 
blur  is  independent  of  the  sampling  and  reconstruction  process  and 
changing  u  merely  causes  the  scene  and  its  image  (Figs,  (la)  and  (lc))  to 
shift  together  along  the  x-axis  without  changing  the  energy  of  their  difference 
Thus  image  blur  is  independent  of  sample-scene  phase.  (This  conclusion  can  be 
verified  mathematically  by  the  substitution  x'  =  x-u  in  Eq.  [4].) 

Sampling  and  reconstruction  blur,  denoted  SR  blur  (Eq.  [5]),  is  just  as 
inevitable  as  image  blur;  however,  unlike  image  blur,  the  amount  and  appearance 
of  SR  blur  depends  upon  the  sample-scene  phase.  Unfortunately,  SR  blur  is 
not  nearly  as  well  understood  as  image  blur,  even  if  the  sample-scene  phase 
effect  is  ignored.  To  think  of  SR  blur  as  just  aliasing  is  to  ignore  the 
effect  of  imperfect  reconstruction;  to  think  of  SR  blur  as  just  interpolation 
(or  resampling)  error  is  to  ignore  the  effect  of  undersampling. 

In  this  paper  the  dependence  of  -t“R  upon  sampling,  reconstruction,  and 
sample-scene  phase  is  analyted.  In  addition,  e^R  is  statistically  character¬ 
ised  in  terms  of  its  mean  and  variance  by  considering  an  ensemble  of  scenes 
f(x-u)  with  u  as  a  random  phase  parameter  uniformly  distributed  between  tero 
and  one.  The  importance  of  this  stochastic  approach  is  that  frequently  in  a 
real  scene  with  many  high  contrast  features,  the  sample-scene  chase  of  each 
is  random  and  equally  likely  to  be  any  number  between  tero  ar.d  one.  Conseauer.t 
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and  therefore  increases  linearly  with  IFOV. 

In  a  discussion  of  image  blur,  s  should  be  interpreted  as  the  scanner 
IFOV.  However,  in  a  discussion  of  SR  blur  it  is  more  appropriate  to  interpret 
s  as  the  number  of  samples  per  IFOV.  Thus,  for  example  s  *  1  corresponds  to 

a  sampling  rate  at  which  consecutive  IFOV's  are  just  contiguous  while  s  =  2 

corresponds  to  a  50%  overlap  in  consecutive  IFOV's.  Most  image  scanning 
systems  operate  with  a  sampling  rate  in  the  range  1  <  s  <  2;  for  mathematical 
simplicity  s  is  restricted  to  this  range  in  all  the  subsequent  spatial 
domain  analysis. 

Sampling  the  ramp  image  (Eq.  [8])  generates  image  samples  g(n-u)  which  are 
all  tero  for  n  <  -1  and  all  one  for  n  >  2.  The  dependence  of  the 
two  remaining  image  samples,  g(-u)  and  g(l-u),  upon  the  parameters  s  and  u 
can  be  determined  from  Eq.  [8].  This  dependence  is  illustrated  in  Fig.  (3). 
For  u  less  than  1  -  0.5s,  g(-u)  =  0.5  -  u/s  and  g(u-I)  *  1;  for  u  between 
1  -  0.5s  and  0.5s,  g(-u)  *  0.5  -  u/s  and  g(l-u)  *  0.3  *  (l-u)/s;  and  for  u 

greater  than  0.5s,  g(-u)  »  0  and  g(i-u)  *  0.5  *  (l-u)./s. 

If  image  reconstruction  is  accomplished  by  linear  interpolation,  the 
associated  interpolation  function  (see  Fig.  (Ie);  is 


r(x) 


t  .  •  ^  Ay 


)1  -  !  x;  <  1 

^0  otherwise 


CIO) 


Figure  (4)  illustrates  the  resulting  reccns true 


d  edge  image  as  well  as  the 


original  edge  and  its  image.  The  shaded  areas  in  Fig.  (4b)  correspond  to  the 
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SR  blur  that  arises  when  1  -  0.5s  <  u  <  0.5s.  Figures  (4a)  and  (4c) 
correspond  to  the  other  two  possibilities,  0  <  u  <  1  -  0.5s  and  0.5s  s  us  1, 
respectively.  It  can  be  shown  that  the  dependence  of  SR  blur  upon  s  and  u  is 


48s 


[(s-2u)2C+2u-s)2  +  (s-2u)2(2-2u-s)2]  0  <  u  <  1  -  0.5s 


SgR(u)  3  \  4gpr  f(5-2u) 2 (2+2u-s)2  +  (s+2u-2)2 (4-2u-s)2]  1  -  0.5s  <  u  <  0.5s 


43s 


■pr  £(s-2u)2 (2-2u+s)2  +  (s+2u-s) 2 (4-2u-s)2J  0.5s  s  u  s  1 


CUD 


As  Eq.  [l 1 ]  indicates,  c2R  is  non-negative  for  all  1  <  s  <  2  and  0  <  u  <  1. 

Figure  (S)  depicts  curves  of  o2R  versus  u  for  various  values  of  s. 

Notice  that  e2R  is  symmetric  about  u  *  0.S  and  that  the  shape  of  these  curves 
depends  dramatically  upon  s;  for  s  =  1  the  curve  is  bowl -shaped  with  maxima 
at  u  *  0,1  and  a  minimum  at  u  =  0.5,  for  s  =  2  the  curve  is  bell-shaped  with 
a  maximum  at  u  =  0.5  and  minima  at  u  ~  0,1.  Equation  [11]  can  be  thought  of 
as  a  transformation,  with  parameter  s,  of  the  uniformly  distributed  random 
variable  u  into  a  new  random  variable  s|R*  The  resulting  probability  density 
function  (pdf)  of  s|R  has  a  very  unusual  shape,  as  illustrated  in  Fig.  (6)  for 
the  case  s  =  1.  The  pdf  is  U-shaped  indicating  that  the  extreme  values  of 
SR  blur  are  significantly  more  likely  than  the  average  value!  The  shape  of 
the  pdf  for  other  values  of  s  is  similar. 

Since  u  is  uniformly  distributed,  the  averages  value  of  SR  blur  is 


e[c1r]  ■ 

and  the  variance  is 

*“[•!*] 


'kiu 


n  ,a 

,  I 


cl.  -  E 


but 


cu 


Cl  2) 


Cl  3) 
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The  variance  is  a  measure  of  the  extent  to  which  depends  upon  u  and  the 


srfect  or 


f  sample-scene  phase  dependence  is  negligible  if  and  only  if  *Var|e|,J 
is  much  less  than  e£c|jJ.  For  example,  from  Fig.  (5),  this  effect  is  negligible 
if  s  »  1.5  but  not  if  s  =  1. 

The  dependence  of  s£c|RJ  and  Var£e|Rl  u?on  s  C3n  determined  algebraical  1) 
from  Eqs.  [ll],  [12]  and  [15].  Specifically,  it  can  be  shown  that 

eBr]  *  1  s  *  5  2  (!J1 

a  surprisingly  simple  result  which  has  been  verified  by  digital  simulation. 

Equation  [14]  illustrates  that  for  an  ensemble  of  edges,  each  blurred  by 

scanning  with  an  ideal  aperture,  sampled,  and  reconstructed  with  linear 

interpolation,  some  SR  blur  is  inevitable;  this  blur  can  be  reduced  by 

increasing  the  sampling  rate,  but  it  cannot  be  eliminated. 

Unlike  Eq.  [14],  the  algebraic  expression  for  Var£egRJ  is  too  cumbersome 

u 

to  have  practical  value  (it  is  an  eighth  degree  polynomial  in  s  divided  by  s  ). 
Instead,  we  present  a  plot  of  SR  blur  variance  versus  sampling  rate,  Fig.  (?) . 
This  variance  changes  by  three  orders  of  magnitude  ever  the  range  i  i  s  <  2 
and  a  pronounced  minimum  occurs  at  a  sampling  rats  of  1.3  samples  per  IrOV. 

It  should  be  emphasized  that  the  previous  analysis  is  valid  only  fer 
1  <  s  $  2.  Intuitively,  the  average  SR  blur  will  continue  to  decrease  in 
the  region  s  >  2  but  not  necessarily  at  the  rate  indicated  in  Eq.  [i-j.  For 
small  sampling  rates,  s  <  1,  it  can  be  shown  that 


Els]  ■  t  -  Ti s  *  7 sI  -  ra 


s  <  1. 


(IS) 
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Fourier  Analysis 

In  order  to  include  the  sample-5cene  phase  effect  into  the  performance 
analysis  of  a  general  image  scanning,  sampling,  and  reconstruction  system  it 


is  necessary  to  formulate  the  problem  in  the  frequency  domain 
13,19 


manner 

oairs 


In  the  usual 

,  we  denote  two  functions  F(x)  and  F(v)  as  Fourier  transform 


F(x)  —  F(v) 


(16a) 


provided 

/•CO 

F(v)  *  |  F(x)e”‘TXV1  dx  (16b) 

J 

and  " 


F(x) 


f  F(v)eJmi  dv 


(16c) 


Since  x  is  normalized  in  units  of  sample  interval,  the  frequency  coordinate  v 
has  units  or  cycles  per  sample  interval.  In  this  coordinate  system  the  Nyquist 
(or  folding)  frequency  is  0.3. 

From  Eqs.  [l],  [2],  and  (3]  it  follows  that 


r(x-u) 


«y  -  2  *uu  i 

" — *  e 


(l'a) 


2( 


..  &  *2"uv i  '  .  ,  ^ 

x-u)  - — -  e  h(v)  f(u) 


:rb) 


and 


&  - 


gT(x;u)  - —  r(v)  £ 


-2^u(v-n)  i 


h(v-n)  f  (v-n) 


(l*c) 


n=-» 


where  f(v)  is  the  (zero  pnase)  scene  frequency  spectrum,  h(v)  is  the  scanner 
OTF,  r(v)  is  the  reconstruction  filter  (i.e.,  the  Fourier  transform  of  the 
interpolation  function),  and  as  before  u  is  the  (random)  sampl e-scer.e  phase 
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parameter.  Equations  [l"a]  and  [l"b]  are  immediate  consequences  of  the  shift 
and  convolution  theorems.  To  derive  Eq.  [17c],  the  conv<"  “ion  theorem  applie 
to  Eq.  [3a]  yields 

gr(x;u)  *—<■  |  ^e'2?rUv  ^  h(v)f(v)J*comb(v)  j  r(v;  ,  (18) 

a  result  which  uses  Eq.  [17b]  and  the  fact  that  the  comb  function  has  the 

18 

curious  property  of  being  its  own  Fourier  transform  .  Convolution  with 
the  comb  function  in  Eq.  [ 1 S ]  produces  the  image  spectrum  replication  in 
Eq.  [17c];  this  replication  is  the  inevitable  result  of  sampling. 

In  the  previous  spatial  domain  analysis  the  scanner  PSF  (Eq.  [7])  and 
interpolation  function  (Eq.  [10])  were  normalised  to  unit  area.  That 
normalization  is  used  throughout  this  paper;  in  the  frequency  domain 
it  means  that  at  zero  frequency 


which  illustrates  (again)  that  image  blur  is  independent  of  sample-scene 
phase.  If  the  scanner  PSF  is  an  impulse,  then  h(v)  =  1  for  all  v  and  image 
blur  is  zero  for  all  scenes.  For  all  practical  PSF's  the  condition  h (%, )  =  1 
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can  not  hold  for  all  v;  at  high  frequencies  h(v)  approaches  cero.  As  Eq.  [20] 
indicates,  the  cause  of  image  blur  is  the  presence  of  significant  scene 

*  A 

energy,  ,'f(v)|2,  at  frequencies  where  h(v)  i  1. 

An  expression  for  SR  blur  analogous  to  Eq.  [20]  can  be  obtained  from 
Eqs.  [5],  [17b]  and  [17c];  however  a  somewhat  different  approach  is  preferable. 
This  approach  is  based  upon  the  observation  that  c|^  is  periodic  in  u  with 
period  one,  and  thus 


eIr<“> 


s  e 
m 


2ru IT  2 


(21a) 


where  the  Fourier  coefficients  are  given  by 


•>  -2Tumi  , 
eSRe  du' 


(21b) 


The  coefficient  a,  has  particular  significance;  it  is  the  average  SR  blur, 
i . e.  ,  from  Eq.  [12] 


[€sr] 


(22) 


In  addition,  it  follows  from  Eqs.  [13]  and  [21a]  that  SR  blur  variance 
is 


J  m=i 


(23) 


In  general  for  m  ±  0  the  Fourier  coefficients  are  complex.  However  because 

e2  is  real,  thev  satisfy  a  =  a*  for  m  =  1,2,3...  where  (•)*  denotes 
sR  -m  m 

complex  conjugate. 

The  explicit  dependence  of  the  Fourier  coefficients  (and  thus  Efe2^]  and 
Var  )  upon  the  scone  spectrum,  scanner  CTF,  and  reconstruction  filter  can 


be  determined  by  first  recognizing  that,  from  Eqs.  [17b]  and  [17c] 


g(x-u)  -  gr(x;u) 


e~2*ruvi  £  -r(v))  h(j-n)  f  (v-n)e2:TUni 

n***  (24) 


where  in  the  summation  6  =  1  if  n  =  0  and  5  =0  otherwise.  A  combination 

n  n 

of  Parseval's  equation,  Eqs.  [S],  [21b]  and  [24]  yields  an  imposing  looking 


expression  for  a 

m 


■n  =  !  21  ^3n-r (u)^ ^o-s-r* (v)^  h(v-n)K'f'j-n)  f  (v-n)  f*(v-n) 


n=-»  n=-» 


e-2^u(ra-n+n) i  du 


L  0 


dv 


(25) 


Fortunately,  because  of  the  integration  with  respect  to  u,  the  only  ncn-tero 
terms  in  the  double  summation  are  those  for  which  n-n  *  m.  With  some 
mathematical  manipulation,  which  is  omitted  for  brevity,  Eq.  [25]  can  be 
reduced  to  its  simplest  form 


m 


r 

i 

—  00 


5m-r(v)-r-(v-m)  +  £  i*(v-n)|2 

]^  =  -ao 


h*  (v ) h  (v  -m)  f  *  (v )  f  (v  -m)  dv 
(26) 


In  summary,  for  a  general  image  scanning,  sampling  and  reconstruction 
system  the  image  blur  is  given  by  Eq.  [20]  and  the  SR  blur  is  given  by  a 
Fourier  series  (Eq.  [21a])  in  the  random  phase  parameter  u  with  coefficients 
given  by  Eq.  [26].  The  mean  and  variance  of  SR  blur  are  given  by  Eqs.  [22] 
and  [23].  All  of  these  quantities  can  be  calculated  by  numerical  quadrature 
provided  values  of  the  scene  frequency  spectrum  f(v),  the  scanner  0 T?  h(v), 
and  the  interpolation  filter  r(v)  are  available. 

As  an  example  of  this  Fourier  analysis,  if  the  scene  is  an  edge  (Eq.  [6 ] ) , 
the  scanner  ?SF  is  an  ideal  aperture  (Eq.  [7]),  and  the  interpolation  function 
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is  linear 

(Eq. 

[10]),  then 

f(u) 

=  4  0  (v)  +  -ry - 

2  i2ttv 

(27a) 

and 

h(v) 

=  sinc(sv) 

(27b) 

r(v) 

=  sinc2(v). 

(27c) 

It  can  be 

shown 

that  the  associated  image  blur  is 

C  tr  — *2 

”1 

!  [1  -  sine (sy) ]2 

J  4t2v2 

(28) 

and  the  SR  blur  is 


CSR(U)  =  E 


[-I*]-  t 

l-  -4  m=l 


B-m  cos(2rrmu) 


where 


'[*]  ■  f[ 


-n)liiS£f 

->  4tt2'. 


1  -  2sinc2(\>)  +  l  sinc^Cv 
n=-® 


and  where  the  Fourier  coefficients  are  real  and  given  by  . 


(su) 


(29a) 


dv  (29b) 


f  -  sine2  (-;)  '  sinc2(v-m)  +  V  sinC*(v-n)| 

L  n=-=  -J 


sinc(S'j)sinc[s  (v-m)]  j, 
4t-v (v-m) 


m  =  1, 


(29c) 


Because  of  the  singularity  in  the  edge  spectrum  (Eq.  [27a])  at  v  =  0,  some  caution 
must  be  exercised  in  the  numerical  integration  of  Eqs.  [23],  [29b],  and  :29c]. 
However,  the  singularities  in  the  integrands  are  removable  and  present  no 
fundamental  difficulties. 
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SR  Blur 

Considerable  insight  into  SR  blur  can  be  obtained  without  resorting  to 
numerical  simulation.  In  particular,  from  Eqs.  [22]  and  [26]  the  average 
SR  blur  is  given  by 


[■*»] 


r 


j 
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e2  (v )  J  h  (v )  f  (v )  [ 2  dv 


where  e2(v)  is  the  non-negative  function  defined  by 

,  ^  30 

e2(v)  =  1  -  2 ,5? [r (v )  ]  -  Y  ;r(v-n)!2 


n=-« 


i 


1  -  r(v)i2  +  £  |  r(v-n)  1 2 

n?0 


(30) 


(31a) 


(31b) 


and  3t  [• ]  denotes  the  real  part  of  a  complex  variable.  Since  e2R  is  a  non¬ 
negative  random  variable,  it  is  identically  oero  if  and  only  if  its  mean  is 
zero.  From  Eq.  [30],  =  0  if  and  only  if  the  integrand  is  identically 

zero  for  ail  \> .  In  other  words,  for  a  general  image  scanning,  sampling  and 
reconstruction  svstem  some  SR  blur  is  inevitable  unless 


e (v )  !h(v)f(v)'  =  0 


iJw 


for  all  v. 

For  a  well  designed  scanning  system  it  is  common  to  use  electronic 
filtering  prior  to  sampling  and  thereby  shape  the  scanner  MTF,  [h(v)j,  so  that 
the  image  spectrum  is  band-limited,  i.e., 

jh(v)  f(v){  =  0  (35) 

for  all  | ')  |  >  v  .  For  such  a  system,  Eq.  [32]  will  be  satisfied  at  high 
frequencies.  However,  below  the  cut-off  frequency,  \>  ,  Eq.  [32]  can  be  satisfied 
only  by  choosing  the  reconstruction  filter  r(v)  so  that  e2(v)  is  zero. 
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What  choice  of  reconstruction  filter  will  satisfy  the  equation 

e-  (v )  =  0  ,  }  v }  <  »;  (34) 

for  al i  | v }  <  \>  ?  Because  of  the  first  term  in  Eq.  [31b],  r(v)  must  be  one  for 
low  frequencies;  because  of  the  remaining  (side-band)  terms  r(u)  must  be  zero 
for  high  frequencies.  In  fact,  in  terms  of  v  ,  the  solution  of  Eq.  [34]  is 
completely  characterized  as  follows.  If  v  >  0.5,  Eq.  [34]  has  no  solution, 
and  if  <  0.5,  the  only  solution  to  Eq.  [34]  is 


r(v ) 


| 1  W\  <  vc 

<  arbitrary  vc  <  'v(  <  0.5 


lo 


f 

V  j 


>  c.s 


CSS) 


This  discussion  establishes  what  amounts  to  a  sampling  theorem"1  for  a 

general  image  scanning,  sampling,  and  reconstruction  system.  Namely,  there 

is  no  SR  blur  (Eq.  [5])  if  and  only  if  all  the  following  are  true: 

(i)  the  image  spectrum  is  band-limited; 

(ii)  the  cut-ofr  frequency  v  is  0.5  cycles  oer  samole  interval  or  less; 

(iii)  the  reconstruction  filter  satisfies  Eq.  [35]. 

Condition  (ii)  represents  the  N'yquist  criteria  for  sufficient  sampling,  i.e.,  it 

determines  an  upper  bound  for  the  sample  interval  size.  Conditions  (i)  and 

(iii)  a.re  idealizations;  a  truly  band-limited  image  spectrum  is  convenient  ts 

assume  but  difficult  to  achieve  and  all  finite  .idth  reconstruction  filters 


are  "non- ideal ,"  i.e.,  they  do  net  satisfy  Eq .  [33].  Thus, 
systems  some  SR  blur  is  inevitable. 

There  are  three  effects  which  contribute  to  SR  blur: 


for  all  practical 


imperfect  recon¬ 


struction  ,  aliasing,  and  sample-scene  phase  dependence, 
of  each,  we  can  combine  Eqs.  [21a]  and  [22]  to  yield 


i 


-4 
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*  £  v5*^  ■ 

From  Eqs.  [30]  and  fclbl,  E  ty  can  be  written  as  the  sum  of  two  terms,  i.e., 

e[£sr]  =  £r  +  £s 


where 


.2  = 


1  -  r  (v  )  |  *  '  h  (v )  f  (v  )  !  2  dv 


(37a) 


(57b) 


pCQ 

s|  =  !  F  y;|r(v-n)!2*]  I  h  (v ) f  (v )  i 2 dv  . 

L  Ln*°  J 


(37c) 


The  term  s2  accounts  for  imperfect  reconstruction;  it  is  analogous  to  image 

blur.  Eg.  [20],  and  it  measures  the  contribution  to  SR  blur  caused  by  the 

presence  of  significant  image  energy,  Jh(v)f(v)|2,  at  frequencies  where  r(v)  4  1 

The  term  £|  accounts  for  aliasing  due  to  undersampling;  it  measures  the 

contribution  to  SR  blur  caused  by  the  presence  of  significant  image  energy 

at  frequencies  where  the  energy  in  the  reconstruction  filter  side-bands, 

^|r(v-n);2,  is  not  cero.  .An  equivalent,  but  more  familiar  expression^4 
n?*0 

for  s2  follows  from  the  identity 


1  H|r(v-n)12  !  h  (v )  f  (v ) '  2 dv 

^  n*0 


2  V 


■  h(v-n)  f  (v-n)  J 2  dv 


Equation  [37d]  demonstrates  the  duality  between  side-band  energy  in  the  recon¬ 
struction  filter  and  image  spectrum.  Both  e2  and  e2  are  independent  of  u 

Kb 

and  correctly  account  for  the  sample-scene  phase-averaged  effects  of  imperfect 
reconstruction  and  aliasing. 

From  Eqs.  [36]  and  [3"a],  e2^  can  be  written 


•2rCu)  =  e2  -  e2  +  b(u) 
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wnere 


>  a  e 

mio  “ 


^TTumi 


'  :  ~L?\' 


(38b ! 


The  term  p  is  the  sample-scene  phase  dependent  contribution  to  SR  blur; 

00 

it  is  a  real  valued  random  variable  with  zero  mean  and  variance  2  7  ia  !2 

m-i  m  ' 

m=i 

Phase  deDendence  is  zero  if  and  only  if  a  =0  for  all  m  =  1,2,...  .  For 
*  m 

a  particular  system  this  phase  dependence  is  best  investigated  by  numerical 
simulation  using  Eqs.  [26]  and  [38b].  However,  one  condition  sufficient  to 
guarantee  that  $  =  0  can  be  obtained  by  inspection;  the  integrand  in  Eq.  [26] 
will  be  zero  for  all  v  if  the  image  spectrum  is  band-limited  with  the  cut-off 
frequency  less  than  0.5.  That  is,  sufficient  sampling  of  a  (truly)  band- 
limited  image  will  eliminate  the  sample-scene  phase  dependence  in  SR  blur. 

Of  course,  some  SR  blur  will  still  be  present  unless  the  reconstruction  filter 
satisfies  Eq.  [55] . 

An  inspection  of  Eq.  [27]  reveals  why  the  example  considered  previously 
in  the  spatial  domain  analysis  exhibited  a  non-zero,  phase  dependent, 

SR  blur  for  all  sampling  rates,  s;  the  image  spectrum  is  not  band-limited 
and  the  reconstruction  filter  is  imperfect.  For  this  example  both 


f  .  .  .]  -> 

s:nc  is; )  - 


1  J  .  .  V5  b.TltlS  .  ]  -  , 

c-  =  j  [l-smc2  (v)  ]-  -tt— -  d-j 


( 59  a) 


r  f  e 


,  J  sincfs';)]2  . 
(v-n)  — — — -  d'J 


(  39b) 


are  positive  for  alls.  In  addition,  it  can  be  verified  numerically  (using 
Eqs.  [23]  and  [29c])  that  VarJt2pJ  >  0  and  thus  the  phase  dependence,  z (u) , 
becomes  snail,  but  does  not  disappear  as  s  becomes  large.  See  also  Fig.  (“) 
which  illustrates  the  magnitude  of  this  phase  dependence  for  1  <  .?  c  2. 
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Reconstruction  Filters 

To  summarize  seme  of  the  previous  discussion,  Eq.  [31b]  defines  e-(v) 

in  terms  of  the  reconstruction  filter  r(v).  As  Eq.  [30]  reveals,  the  cause 

of  SR  blur  is  the  presence  of  significant  image  energy  at  frequencies  where 

e2(v)  is  not  zero.  Just  as  the  scanner  OTF  can  be  interpreted  as  a  pre- 

sampling  filter  acting  upon  the  scene  Wiener  spectrum,  jf(v)[2,  to  produce 

image  blur  (Eq.  [20]],  the  reconstruction  filter  can  be  interpreted  as  a 

post-sampling  filter  acting  upon  the  image  Wiener  spectrum,  jh(\>)  f(v) !  2  ,  to 

produce  average  SR  blur  (Eq.  [30]).  In  other  words,  ]  1  -h ('j ')  { 2  is  to  image 

blur  as  e2(v)  is  to  average  SR  blur.  It  should  be  emphasized  that  e2(\>) 

is  not  just  jl-r(v)l2;  sampling  causes  the  appearance  of  the  side  band  terms 

2[r(v-n)[2  in  Eq.  [3ib].  If  r(x)  is  the  ideal  interpolation  function  sinc(x) , 
n^O  „ 

then  r(\>)  *  rect(v)  and  e2  (■•;)  =  2  -  2rect(v).  What  does  e2(\>)  look  like  for 
some  common  (non-ideal)  reconstruction  filters? 

Table  I  lists  the  characteristics  of  three  common  digital  interpolation 
functions,  nearest  neighbor,  linear,  and  cubic  and  the  ideal  interpolation 
function,  sine.  There  are  actually  several  variants  of  the  cubic  interpolator 
in  the  literature,  each  derived  as  a  finite-width  approximation  to  sinc(x); 
the  form  we  have  chosen  is  due  to  3emsrsin‘"’ .  Figure  (3)  is  a  plot  cf  e2(-;} 
versus  v  for  the  four  interpolation  functions.  A  comparison  of  the  relative 
magnitude  of  the  four  curves  in  the  region  0  s  v  s  0.3  yields  the  conclusion 
that  for  any  band-limited  and  sufficiently  sampled  image:  (i)  the  ideal 
interpolator  causes  no  SR  blur;  (ii)  of  the  three  common  digital  interpolation 
functions,  cubic  causes  the  least  SR  blur  and;  (iii)  nearest  neighbor  inter¬ 
polation  causes  the  most  SR  blur.  (Unfortunately,  it  is  necessary 
to  qualify  conclusion  (ii)  because  cf  the  small  hump  on  the  cubic  curve  at  low 
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frequency.  It  would  be  mathematically  possible  to  construct  a  band-limited 
and  sufficiently  sampled  image,  with  virtually  ail  its  energy  at  ■,  =  0.0S, 
for  which  linear  interpolation  causes  slightly  less  SR  blur  than  cubic.) 

Since  the  e2(v)  curves  cross  above  the  Nyquist  frequency  (0.S),  if  the  image 
is  not  band-limited  and  sufficiently  sampled,  conclusions  fi),  (ii),  and  (iii) 
are  not  necessarily  valid.  In  fact,  in  this  case  conclusion  (i)  is  invalid 
and  the  ideal  interpolator  may  cause  more  average  SR  blur  than  some  of  the 
others! 

The  example  we  have  considered  throughout  this  paper--eage  reconstruction-- 
has  an  image  spectrum  which  is  not  band-limited.  For  this  example  the  average 
SR  blur  is 


-  r  ->  n  r  e2  Cv )  sine2  (s\>)  , 

E  LeSRJ  "  j  dv 


(40) 


where  e2(v)  is  given  by  Eq.  [51]  and  the  terms  r(v)  and  £  !r(v-n)|2  are  listed 
in  Table  I  for  each  interpolation  function.  Equation  [40]  has  been  evaluated 
numerically  for  various  values  of  the  sampling  rate,  s,  and  the  results  are 
illustrated  in  Fig.  (9).  Of  the  three  common  digital  interpolators,  nearest 
neighbor  is  clearly  the  worst  choice  (as  expected)  and  cubic  is  the  best 
choice  except  at  low  sampling  rates,  s  %  0.6.  In  this  example,  when  s  is 
less  than  0.3,  the  image  spectrum  has  so  much  energy  above  the  Nyquist  frequency 
that  the  ideal  interpolator  actually  causes  more  average  SR  blur  than  either 
linear  or  cubic.  In  fact,  cubic  causes  less  average  SR  blur  than  sine  for 
values  of  s  as  large  as  1.1. 

Table  I  lists,  for  each  interpolator,  a  closed  form  expression  for  the 

infinite  series  2_.  ' r(v-n)|2.  This  expression  car.  be  used  to  simplify 

n*-® 

considerably  the  numerical  evaluation  of  e2(v)  and  is  based  upon  the  identity 
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33  30 

£  '  rtv-n)'2  =  (r*rl(0)  *  2  ^  (r*r)  (n)  cos(2-rnv)  (41) 

n»-»  n = 1 

Equation  [41]  is  valid  for  any  interpolation  function  provided  it  is  even, 
i.e.,  r(x)  *  r(-x)  for  all  x.  For  all  practical  interpolation  functions, 

(r*r)(n)  will  be  tero  except  at  small  n  and  the  infinite  series  in  Eq.  [4 1 ] 
can  be  evaluated  in  closed  form.  The  validity  of  Eq.  [41]  follows  from  the 

observations  that  if  r(x)  is  even,  r(v)  will  be  real,  r(x)*r(x)  will  be  even, 
and 

X 

23  |r(v-n)|2  =  r2(v)*  comb(v)  =  21  (r*r)(n)e  2i:nvi_  C42 ) 

n=-<»  n=-« 

The  last  equality  in  Eq.  [42 J  can  be  verified  by  taking  the  Fourier  transform 
of  each  side.  Ir.  Eqs.  [41]  and  [42],  (r*r)  (n)  represents  r(x)*r(x)  evaluated 
at  x  =  n. 

The  assumption  that  the  interpolation  function  r(x)  is  even  is  not  restricti 

it  is  difficult  to  imagine  any  situation  in  which  r(x)  would  be  constructed 

otherwise.  It  can  be  shown  that  if  r(x)  is  even,  then  r(;),  23! *(v-n) i 2 » 

n 

and  e2(v)  will  all  be  even  functions  of  v. 


Discussion 


The  emphasis  throughout  this  paper  has  been  on  an  analysis  of  the  degradation 
e|^(u),  caused  by  image  sampling  and  reconstruction.  The  analysis  is  unique  in 
that  it  includes  the  sample-scene  phase  as  an  explicit,  random  parameter. 

Results  established  in  the  frequency  domain  provide  a  complete  characterisation 
of  c~^(u)  as  the  sum  of  two  terms,  u|e|^J  and  i(u).  The  first  term  accounts 
for  the  mean  effect  of  sampling  and  reconstruction  averaged  over  all 
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sample-scene  phases,  the  second  term  accounts  for  variations  about  this  mean. 

If  the  image  is  band-limited  and  sufficiently  sampled,  the  variations  about  the 
mean  are  suppressed,  i.e.,  c (u)  =  0.  If,  in  addition,  the  reconstruction 
filter  is  ideal,  the  mean  is  rero  and  SR  blur  disappears. 

Table  II  is  a  logic  flow  diagram  which  summarizes  the  results  of  the 
frequency  domain  analysis.  To  avoid  unnecessary  clutter  in  this  diagram,  only 
a  band-limited  image  is  considered  and  the  notation 

g(v)  =  h(v)f(v)  (45a) 


v(v  ;m) 


—  -  A  A  A 

;  r (v-n)  j  2  -  r(v)  -  r*(v-m) 

n=-<*> 


(43b) 


is  introduced  to  simplify  some  equations.  Since  all  spatial  and  frequency 
variables  in  this  paper  are  referenced  to  a  coordinate  system  in  which  the 
sample  interval  has  unit  length,  all  frequencies  (\>)  in  Table  II  have  units  of 
cycles  per  sample  interval  and  the  N'yquist  frequency  is  0.5. 

To  better  understand  the  results  of  this  paper,  it  is  important  to 
appreciate  the  inverse  relation  between  the  parameters  v  and  s.  For  a  general 
scanner  PSF  (not  just  an  ideal  aperture,  Eq.  [7]),  s  represents  the  width  of 
the  IFOY  relative  to  the  sample  interval  or,  equivalently,  the  number  of 
samples  per  IFOV.  Just  as  s  is  proportional  to  the  effective  width  of  the 
scanner  PSF,  the  cut-off  frequency,  v  ,  is  proportional  to  the  effective 
width  of  the  scanner  CTF.  Because  of  the  scaling  property  of  Fourier  transform 
pairs,  it  follows  that 


s  v. 


constant. 


(44) 
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Thus,  for  example,  reducing  the  cut-off  frequency  by  SO -a  is  equivalent  to  doubling 
the  sampling  rate. 

From  Eq.  [30]  it  can  be  seen  that,  for  a  fixed  sample  interval  size,  there 

are  two  ways  to  reduce  SR  blur.  One  way  is  to  shape  the  reconstruction  filter 

so  that  e2  (v)  is  small  at  those  frequencies  where  [h(v)f(v)!~  is  large.  As 

Fig.  (S)  indicates,  this  technique  of  post-sample  filtering  can  be  quite  effective, 

but  only  if  the  IFOV  width,  s,  is  large  enough  so  that  v  <  0.5.  Reducing  SR 

blur  by  shaping  the  reconstruction  filter  does  net  effect  image  blur.  The  other 

* 

way  to  reduce  SR  blur  is  to  shape  the  scanner  MTF  so  that  jh(v)!:  is  small  at 
those  frequencies  where  e2  (v )  |  f  (v )  I  -  is  large.  This  technique  of  pre-sample 
filtering  definitely  does  effect  image  blur.  In  fact,  there  is  a  trade-off: 
reducing  SR  blur  by  shaping  the  scanner  MTF  inevitably  increases  image  blur. 

Figure  (9}  illustrates  the  trade-off  between  image  blur  and  SR  blur  in  the 
special  case  of  an  edge  scanned  with  an  ideal  aperture.  For  a  fixed  sample 
interval  site,  as  s  increases  the  IFOV  is  broadened,  the  effective  cut-off 
frequency  is  reduced,  and  high  frequencies  are  suppressed.  The  result  is  an 
increase  in  image  blur  and  a  decrease  in  SR  blur. 

Finally,  it  is  important  to  note  that  the  total  blur  (i.e.,  the  square  of 
the  total  radiometric  error)  associated  with  the  entire  process  of  scanning, 
sampling  and  reconstruction  is 


*•30 

s2(u)  =  j  [f(x-u)  -  g„(x;u)J2dx  , 

“  CO 

This  sample-scene  phase  dependent  blur  is  not  just  the  sum  of 
However,  from  the  triangle  inequality,  it  follows  that 


(43) 


and 


e2(u) 


‘IrW 


(46a) 

(46b) 
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Therefore,  the  results  of  this  paper  can  be  used  directly  to  establish  an  upper 
bound  on  -:2  and  E[e2]. 


(a)  SCENE:  f(x-u) 
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Table  II.  The  calculation  of  SR  blur.  The  image  spectrum,  g(\>) ,  is 
assumed  to  be  cero  for  all  [v[  >  vc. 


INPUT 

SCENE  SPECTRUM  f(v) 

SCANNER  OTF  h(v) 

CUT-OFF  FREQUENCY  v  (CYC/SAMPLE  INTERVAL) 
RECONSTRUCTION  FILTER  f(v) 


CALCULATE  BAND  LIMITED  IMAGE  lh(v)?(v)  J v j <v c 

SPECTRUM  g(v)  ^0  otherwise 

e2(y)  =  l-25?[r  (v)]+T|  r(v-n)  j2 

n 

v(v;m)  =  £j  r(v-n)  |  2- r  (v) -r*  (v-m) 

_ _ n  . . . - 

esr(u)  "  E[eiRh(u) 


vc  <  °.5 

(SUFFICIENT 
SAMPLING)  . 


E{s|rJ  =  fCe2(v)|g(v)|2dv 


i  \  r  2*umi] 

<?(u)  -  2  i  ^[ame  j 

m=l 


v(v ;m) g« (v) g(v-m) dv 


/'DOESX 

r(v)  \ 
SATISFY  \ 
.Eq.  (351  (IDEAL 

\  F,L?T£V 


rf,2  J 

i'SR! 


r2  =  e  l’c.2  1 
'SR  “i'SR! 

r>c 

=  e-  {-j  )  g  (')  i  ■  2dv 


($(u)  =*  0) 


c-  =  0 

CSR  J 
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LIST  OF  FIGURE  CAPTIONS 


Fig.  1.  The  image  scanning,  sampling,  and  reconstruction  process. 

Fig.  2.  (a)  The  blurring  introduced  by  scanning.  (b)  The  additional 

blurring  introduced  by  sampling  and  reconstruction. 

Fig.  5.  Ramp  image  sample  dependence  upon  the  sample-scene  phase  parameter. 

Fig.  4.  SR  blur  dependence  upon  the  sample-scene  phase  parameter.  Three 
distinct  situations  are  depicted  corresponding  to  (a)  0<u<l-0.5s, 
(b)  1-0. Ss<:u<0. 5s,  and  (c)  0.5<u<l,  respectively. 

Fig.  5.  The  dependence  of  SR  blur  upon  the  sample-scene  phase  parameter. 


The  curves  correspond  to  s=1.0,  1.25,  1.50,  1.75,  and  2.00,  respectively. 


For  s=1.50,  c2R  is  nearly  constant. 

Fig.  6.  The  probability  density  (vertical  axis)  of  for  the  case  s=l. 

The  extreme  values  of  e2R  are  significantly  more  likely  than  the 
average  value,  E [ sSR  I ' 

Fig.  7.  SR  blur  variance  versus  sampling  rate  (samples/'IFOV)  for  edge  recon¬ 


struction. 

Fig.  S.  The  function  e2 (v)  for  three  common  interpolators,  nearest  neighbor, 

linear,  cubic  and  the  ideal  interpolator,  sine.  The  insert  illustrates 
the  small  hump  at  low  frequencies  for  cubic. 

Fig.  9.  Average  SR  blur  as  a  function  of  sampling  rate  (samples/ I FOV)  for  the 


interpolators,  nearest  neighbor,  linear,  cubic  and  sine.  The 
image  is  an  edge  scanned  with  an  ideal  aperture. 
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LIST  OF  TABLES 

Table  I.  Characteristics  of  three  common  digital  reconstruction  functions, 

nearest  neighbor,  linear,  cubic,  and  the  ideal  interpolator,  sine. 

Table  II.  The  calculation  of  SR  blur.  The  image  spectrum,  g(v),  is  assumed  to 

be  zero  for  all  ' u !  >  v  . 
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Research  Objectives 

Interframe  data  compression  with  optical  computations  has 
been  demonstrated  under  previous  support  of  AFOSR.  Interframe 
data  compression  has  proven  more  elusive,  principally  because  of 
the  complex  logic  operations  required  in  interframe  comparisons. 

The  objective  of  this  research  is  to  demonstrate  that  architectures 
for  interframe  compression  can  be  conceived  which  combine  computa¬ 
tions  appropriate  with  digital  processors. 

Appendix  Contents 

The  contents  of  this  appendix  consists  of  a  review  of  the 
general  problem  of  interframe  compression  plus  a  proposed 
architecture  for  an  opti cal/digi  tal  system  for  same,  plus  a  paper 
demonstrating  the  results  achieved  in  simulations  of  one  such 
method. 


I.  INTRODUCTION 


A  hybrid  optical/digital  system  architecture  has  been  shown  in  previous 
research  to  be  a  feasible  and  viable  approach  to  interframe  coder  design 
[1,2].  One  of  the  definite  advantages  of  optical  computation  over  digital 
computation  is  the  speed  and  ease  of  frame  differencing  in  terns  of  the 
inherent  image-piane-to-image-plane  parallel  processing.  The  temporally 
adaptive  hybrid  optical/digital  system  proposed  in  Ref.  [2]  was  designed  to 
use  the  above  advantage,  as  well  as  the  temporal  redundancy  between 
successive  image  sequences.  However,  it  is  not  quite  spatially  adaptive  due 
to  the  fact  that  only  the  interpolated  differential  pulse  code  modulation 
(1DPCM)  spatial  compression  system  has  been  used  where  large  movement 
greater  than  the  predetermined  threshold  is  involved  between  image 
sequences.  Therefore,  the  above  scheme  can  be  greatly  improved  by 
considering  the  following  factors: 

(1)  An  adaptively-coordinated  spatial/ temporal  subsampling  and 
interpolation  between  the  spatial  and  temporal  compression 
subsystems  are  feasible  by  using  programmable  resolution 
provided  by  focal  plane  detector  technology. 

(2)  Motion-detection/estination  schemes  could  be  integrated  into 
a  hybrid  optical/digital  architecture. 

(3)  More  reliable  movement  activity  measures  (image  derivative 
energy)  [3]  could  replace  the  normalized  mean  square  error 
(NI1SE)  of  the  frame  difference,  which  could  be  used  for  the 
decision  mechanism  of  an  adaptive  compression  system. 


(A)  A  highiy  parallel  digital  processing  scheme  using  multiple 
microprocessors  could  be  integrated  into  the  hybrid 
optical/digital  architecture  in  order  to  offset  the  slow 
processing  speed  of  sequential  digital  subsystems  as  opposed 
to  the  high-speed  optical  subsystem.  Also,  the  problem  of 
efficient  interfacing  between  the  optical  and  digital 
subsystems  should  not  be  overlooked. 

This  paper  addresses  the  above  suggestions  for  the  improvement  of  the 
hybrid  optical/digital  approach  for  interfraue  image  data  compression  by 
introducing  the  adaptive  decision  mechanism  based  on  image  derivative  energy. 
In  Chapter  II,  the  state  of  the  art  of  the  adaptive  digital  interframe 

predictive  coding  is  reviewed.  In  Chapter  III,  we  propose  the  use  of  image 

> 

activity  measurement  for  adaptive  multimode  coding,  which  is  designed  to 
integrate  high-speed  optical  processing  for  the  computation  of  spatial  and 
temporal  gradients  and  energies,  and  flexible  digital  processing  for 
controlling  multimode  coding.  In  addition  to  the  system  configuration  and 
algorithmic  explanation,  the  current  technological  trends  are  briefly 
described  to  justify  parallel  processing  and  hybrid  optical/digital 
processing. 

Finally,  the  striking  similarity  between  the  parallel  image  processor  in 
the  proposed  adaptive  multimode  coder  and  the  parallel  processor  ?ASM 
(partitionable  SIMD/MIMD  system)  has  been  pointed  out.  Thus,  it  is  inferred 
that  the  parallel  hybrid  optical/ digital  processing  approach  provides  a  much 
broader  conceptual  framework  fer  computer  vision  systems.  Detailed 
architectural  consideration  will  be  left  for  further  research. 
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II.  ADAPTIVE  DIGITAL  INTERFRAME  PREDICTIVE  CODING:  A  REVIEW 

There  have  been  numerous  excellent  surveys  concerning  digital  picture 
coding  14-91.  Rather  than  repeating  the  above  surveys,  we  selectively 
consider  the  recent  development  related  to  adaptive  interframe  predictive 
coding. 

There  are  various  adaptive  strategies  that  have  been  developed  for 
spatial  compression  algorithms,  such  as  differential  pulse  code  modulation 
(DPCM),  etc.  The  typical  examples  are  adaptive  prediction  and  quantitation 
approaches.  The  former  is  based  on  the  nonstationary  statistical  image 
model  and  assumption,  which  are  more  reliable  in  characterizing  the  edges  or 
outline  of  the  real  world  images  than  the  first  order  Markov  image  model. 
Thus,  the  prediction  error  at  the  edge  between  the  adaptive  predictor  and 
the  original  pixel  is  made  smaller  than  the  nor.adaptive  prediction  error. 
The  latter  takes  advantage  of  the  insensitivity  of  the  human  eye  to  the 
quantization  error.  That  is,  the  quantization  steps  and  levels  are 
adaptively  selected  in  such  a  way  that  the  introduced  quantization 
distortion  is  not  perceptable  to  the  human  visual  system. 

Interframe  compression  algorithms  also  take  advantage  of  the  temporal 
characteristics  of  the  human  visual  system,  in  addition  to  the  above  spatial 
characteristics  of  human  visual  perception  and  the  statistical  redundancy 
represented  by  an  image  source  model.  A  recent  study  [9j  on  the  spatio- 
temporal  response  of  the  human  visual  system  revealed  that,  at  high 
temporal  frequencies,  spatial  contrast  sensitivity  is  reduced.  This  indicates 
that  the  sensitivity  of  the  human  visual  system  is  low  to  fine  spatial 
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detail  in  rapidly-moving  objects,  and  at  high  spatial  frequencies  there  is  an 
overall  decrease  in  flicker  sensitivity.  For  the  nontracking  task,  which  is 
typical  in  teleconference  or  picture  phone  environments,  the  human  visual 
system  can  tolerate  a  loss  of  spatial  resolution  in  the  moving  area 
(reported  by  Miyahara  [10)).  These  observations  strongly  suggest  that  the 
interframe  coder  can  classify  the  parts  of  the  spatial  source  image 
sequences  and  adaptively  apply  different  efficient  coding  schemes  to  them, 
in  accordance  with  the  different  spatial  and  temporal  resolutions  required 
by  the  human  viewer  for  the  appropriate  visual  perception.  The  simplest 
coding  schemes  are  the  frame  replenishment  method  using  frame  repeating¬ 
dropping,  spatial  and  temporal  subsampling  and  interpolation. 

The  more  sophisticated  methods  consist  of  extracting  the  actual  motion 
displacement  between  frames.  This  motion  estimation  approach  is  still  in 
its  infancy  and  most  of  the  past  work  has  been  restricted  to  two- 
dimensional  motion,  especially  translation  [9,11).  There  have  been  four 
major  approaches  to  estimating  tw'o-dimer.sicnal  translation:  (1)  the  Fourier 
method,  (2)  the  method  of  differencing,  (3)  the  temporal-spatial  gradient 
technique,  and  (4)  matching.  In  the  context  of  TV  signal  transmission,  the 
latter  three  have  been  selected  for  real-time  implementation.  The  basis  of 
movement-compensated  predictive  coding  is  that,  if  the  displacement  field  of 
the  moving-image  sequence  is  known,  then  a  very  good  prediction  of  the 
current  field  in  the  moving  image  sequence  can  be  obtained  by  shifting  and 
interpolating  the  parts  of  the  previous  field  (or  frame)  which  have  moved 
[9].  This  general  approach  is  also  simplified  into  two  methods  for  real-time 


impie  mentation. 


Cr.e  method  is  the  Pel-block  displacement 


66 


estimation/corapensation  technique  [12-18]  and  the  ether  is  the  Pel-recursive 
displacement  estimation/compensation  technique  [19,20,25]. 

In  the  former,  the  scene  is  divided  into  rectangular  blocks,  and  a 
single  shift  is  estimated  for  each  block.  The  block  may  or  may  not  be 
segmented  into  changed  and  stationary  areas.  The  previous  frame  (or  field) 
is  shifted,  and  interpolated  if  nonintegral  displacement  estimates  are 
allowed,  to  form  the  prediction  for  that  block  [9]. 

The  latter  technique  updates  the  displacement  estimate  (which  may  or 
may  not  be  obtained  by  Pel-block  motion  estimation)  at  each  pixel.  This 
update  is  based  only  on  previously  transmitted  pels,  so  that  no  explicit 
displacement  estimate  need  be  transmitted. 

There  are  two  methods  for  estimating  translational  displacements:  (1) 
the  correlation  or  matching  technique  [  1.2,13,18],  and  (2)  the  temporal-spatial 
gradient  technique  [10,14,17,21,22,23].  Given  an  object  in  translational 
notion  with  velocity  VKV^,!^),  the  image  luminance  satisfies 

U(X,t)  =  U[X-V(t-t0J,t0]  (II.l) 

for  an  arbitrary  reference  time  tg.  In  particular,  if  t  is  the  time  interval 
between  frames 


yjc.tg+T) 


>c0.’ 


(11.2) 


where  d=TV  is  the  displacement  which  occurs  in  one  frame  interval. 
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II.  1.  The  Correlation/Matching  Technique 

The  basic  idea  of  this  method  is  to  choose  d_  from  a  set  of  possible 
displacements  3  to  minimize  some  measure  of  the  difference  between  the 
current  frame  and  the  previous  frame  displaced  by  d.  Thus,  d_  is  the  value 
which  solves 


min  l  N[u(x,t0+7j-u(x^d,t0j]  (II.3) 

dc3  XaMA 

where  MA  denotes  the  moving  area  and  N(*)  is  a  nonnegative  increasing 
function  which  serves  as  a  distance  measure. 

£i^UHtlOn  (ll*  3)  indicates  that  a  large  computation  is  required,  since  it 
should  be  evaluated  for  each  dz3  .  Thus,  there  is  a  trade-off  between 
accuracy  and  computational  load.  In  practice,  only  a  relatively  small 
number  of  integral  displacements  have  been  used. 


II. 2.  The  Temporal-Spatial  Gradient  Technique  [14,21,22,23] 

The  standard  procedure  to  minimize  the  difference  between  l:(X,t)  and 
U(X,t)*U[X^y(t-t0j,to]  is  to  take  the  derivative  of  the  squared  error 


z=[U(X,t)-U(X  ,t)]~  and  equate  it  to  zero, 
between  the  temporal  derivative  and  the 
object  as  shown  in  Ref.  [24J: 


Since  there  exists 
spatial  gradient  of 


u(-:. 


), 


a  relation 
the  moving 


(II. 4) 
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approximating  the  derivatives  by  finite  differences,  we  have 

AT(JC,t)  =  -  V'Ax(X,t). 

Thus,  the  velocity  estimate  is  obtained  by  choosing  V,  which  giv 
mean-square  fit  to  (II.5)  over  the  moving  area,  i.e.,  which  minimi: 

C  =  l  i  AT(X,t)  +  V*Ax(X,t)j2. 

XeilA  ~ 

By  setting  AyC=0,  we  obtain 

V  -  -I  l  uxCX,t)AxCX,t)Tr1  _  uT(X,t)Ax(X,t). 
[XeMA  ~  j  •  XeMA  ~ 


Limb  and  Murphy  [10,17]  have  described  a  simplified  version  of 
estimator  as  follows: 


V  = 


1  Ax(_X,t)  sign]  Av.  (X,t)j/  2. 

XeMA  '  /  XeMA 


XeMA 


L  Ax(X,t;  sign!  Ax  (_X.t)]/  \ 

f  XeMA 


Ax2(_x,t)! 


whe-e  X=*vXj,X-j). 


CII.5) 

es  the  best 
:es 

CII.6) 


(II.7) 

the  above 


(11.8) 


A  survey  of  these  two  methods  is  given  in  Ref.  111!- 


The  gradient 


methods  in  general  require  fewer  computations  than  the  matching/correlation 
technique  for  a  given  level  of  accuracy.  This  is  particularly  critical  for 
real-time  application  in  image  sequence  (interframe)  coding.  However,  the 
gradient  method  is  inaccurate  for  large  displacements  (above  several  pixels 
per  frame)  where  there  is  a  significant  amount  of  detail  due  to  the  fact 
that  the  temporal  derivative  in  Eq.  (II.4)  is  approximated  by  the  frame 
difference.  This  problem  can  be  alleviated  by  introducing  a  temporal 
recursion.  In  ether  words,  the  corresponding  estimate  for  the  previous 
frame  is  updated  using  the  same  algorithm  with  the  Pel-block  displacement 
estimation  technique,  but  the  frame  difference  is  replaced  by  the  displaced 
frame  difference.  The  resulting  estimator,  called  the  directional  gradient 
method  or  the  Pel-recursive  displacement  estimation  technique,  is  accurate 
and  robust  to  noise  for  large  area  translation,  but  it  cannot  respond  to 
rapid  spatial  changes  in  the  displacement  field. 

Netravali  and  Robbins  [19]  have  proposed  a  coder  in  which  the 
displacement  estimate  is  recursively  updated  at  a  spatially-neighboring 
point  as  follows:  Given  the  displaced  frame  difference,  with  displacement  d_, 

D(X,t,d)  =  U(JC,t)  -  U(X-_d,t-T),  (II.9) 

and  if  the  object  has  translated  by  an  amount  d_  in  one  frame,  then 
D(X,t,d)=0.  The  basis  of  the  pel-recursive  estimator  is  to  adjust  d_  at  each 
pel  in  order  to  reduce  |DOC,t,d)j.  This  is  done  by  changing  d_  in  the 
direction  of  the  negative  gradient,  that  is. 
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ii_1  -  |  7d[D(X,ttdi-lJ]< 


(n.io) 


where  d_*  is  Che  displacement  estimate  at  pel  i  of  some  scanning  sequence. 
Inserting  Eq.  (.11.9)  into  Eq.  (II. 10)  yields 


-d1'1 


eD[X,t,di-ljVxU^i-1,t-Tj. 


(11.11) 


Since  U(X,t)  is  defined  only  at  the  sampling  point,  interpolation  must  be 

A  A 

used  to  evaluate  D(X,t,_d)  and  VxU(X^d  ,t-t).  They  found  it  adequate  to  use 
linear  interpolation  to  evaluate  U(X^d  ,t-x),  and  7xU'X-[ dj ,t-x J  to  estimate 

A  A 

the  gradient,  where  (dj  is  d  rounded  to  the  nearest  grid  point.  This 
algorithm,  where  £“1/1024,  further  reduced  the  bit  rate  by  1  to  1.5  bits  per 
pixel  for  the  dispiaced-f rame  difference  with  respect  to  simple  frame 
difference.  However,  the  Pel-recursive  algorithm  has  limited  spatial 
convergence  to  ensure  stability  and  is  more  sensitive  to  noise. 

The  Pel-recursive  algorithm  is  also  accompanied  by  the  problem  of 
addressing  the  moving  and  stationary  areas  and  coding  the  movement- 
compensated  prediction  error  in  moving  areas.  Upon  successful  movement 
compensation,  most  of  the  quantized  prediction  errors  will  be  zero,  thus,  use 
of  efficient  run-length  coding  is  justified.  Also,  since  there  are  three 
states  (stationary  area,  insignificant  compensated  prediction  errors,  and 
significant  compensated  prediction  errors),  an  identification  for  the  state 
of  the  next  run  length  must  be  provided  in  addition  to  the  run  lengths.  An 


alternative  to  this  three-state  encoding  is  to  run-length  code  only 
significant  and  insignificant  prediction  errors.  The  decision  as  to  which 
prediction  to  use  (i.e.,  previous  frame  or  displaced  previous  frame)  can  be 
made  by  using  previously  transmitted  information,  so  that  no  over  head 
information  need  be  transmitted. 

The  Pel-recursive  displacement  estimation  has  been  extended  to  the 
transform  domain  in  Ref.  [25J.  The  algorithm  works  recursively  in  such  a 
way  that  the  displacement  estimates  are  updated  from  two-dimensional 
unitary  transform  coefficients  to  two-dimensional  unitary  transform 
coefficients.  Each  image  is  partitioned  into  blocks  of  size  N'“N'r  *N0,  to 
obtain  the  N  transform  coefficients  for  each  qth  block  subimage, 

Cn(q)  =  UT(Xq,tjPn 

Xq  -  [Xlq,X2qjT  (11.12) 

where 

Xq  =  iXiq,X2q]”  denote  the  coordinates  of  the  upper  left-hand  Pel 
of  the  qth  block 

U(Xq,tj  =  the  column  scanned  vector  of  intensities  of  block  q 
Cn(q)  =  the  nth  transform  coefficient 

=  the  nth  linear  transform  basis  vector. 
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is  the  nth  sufficient  in  the  previous  displaced  frame.  When  Cn( q)  is 
predicted  by  Cn(q,_D),  the  prediction  error  is 

C(q,D)  =  Cn(q)  -  C„Cq,D).  (11.14) 

The  squared  prediction  errors  are  minimized  by  a  steepest  descent  method. 
The  resulting  recursive  formula  is 

D.n+l<<l>  =  -  I  '7_Dn(q)eh(q.Dn<<l)j* 

Jones  and  Rashid  [20]  have  proposed  residual-recursive  displacement 
estimation,  which  utilizes  adaptive  estimation  techniques  to  vary  certain 
parameters  of  adaptive  hybrid  picture  coding  (AHPC)  as  the  statistics  of  the 
image  change,  in  AHPC,  a  one-dimensionally  transformed  S(k)  in^a  raster  scan 
fashion  is  predicted  by 

? 

SGO  =  :  a^(k)S(k-i)  +  e(k),  (11.15) 

i-1 

where  the  predictor  coefficient  a^(k;  is  a  piece-wise  constant.  If  a^(k)  is 
estimated,  the  residual  sequence  elements  are  giver,  by 

P  „ 

e(k)  =  5(k)  -  .■  ai(k)S(k-i),  (11.16) 

i=l 

where  e(k)  is  the  difference  between  the  prediction  of  the  transformed 
picture  element  and  the  actual  value.  Assume  that  the  image  is  divided  into 


blocks  of  Nr  rows  by  Nc  columns.  Then  each  residual  in  che  present  frame 


may  be  represented  by 


^n  ^  ^  J  ~  ^  ^  j  **n  "  ^  (jlc  > b )  3-ji  * 


(n.17) 


where  '^(Xq.t)  ;n  is  a  linear  prediction. 

Similarly,  the  displaced  previous  frame  value  for  this  residual  is 


en(q,_D)  “  i*'T fXq -JD  ,t-T'i4n  -  UT(Xq-D,t-T )  ©n. 


-  iiT/v  _r 


(11.16) 


A  displaced  residual  difference  may  be  formed  by 


DRDn(q,D)  =  en(q)  -  en(q,D). 


(11.19) 


The  algorithm  attempts  to  minimize  the  squared  displaced  residual  difference 
in  a  residual  recursive  steepest  descent  procedure  as  follows: 


En+1<^  =  -  j  rVn(q)j_DRDn(q,D)j2. 


(H.20) 


Tlris  operation  may  be  simplified  by 


Dn+i(q)  =  Dn(q)  -  cDRD(q,D(q))en(Xq-_D  ,T~t)).  (H.21) 


They  have  shown  that  this  algorithm  is  more  stable  in  the  region  e  =» 
lU-4  than  the  coefficient  recursive  displacement  estimation  algorithm. 


L_ 
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ill.  A  SPATIO-TEMPCRALLY  ADAPTIVE  MULTIMODE  CODING  UNDER  PARALLEL 
HYBRID  OPTICAL/ DIGITAL  PROCESSING  CONCEPT 

Ill.i.  Adaptive  Multimode  Coding  Using  Image  Activity  Measure 

As  discussed  :'r.  the  previous  chapter,  the  state-of-the-art  interframe 
coding  algorithms  have  become  very  sophisticated  at  the  cost  of 
implementation  and  system  complexity.  However,  the  compressed  bit  rates 
remain  beyond  the  capabilities  of  conventional  digital  transmission  channels, 
such  as  telephone  links.  The  common  characteristic  of  most  interframe 
coding  implementation  is  the  need  for  a  frame  buffer  for  storing  and 
smoothing  the  image  information  to  be  transmitted.  Since  the  size  of  the 
buffer  is  limited  by  constraints  relating  to  cost  and  perceptual  delay 
effects,  peaks  of  activity  in  the  image  signal  can  cause  this  buffer  to 
overflow.  To  avoid  this,  it  is  necessary  to  adapt  the  coder  operation  to  the 
amount  of  motion,  that  is,  to  deliberately  degrade  the  image  quality  in  a 
gradual  and  graceful  manner  as  the  image  activity  increases.  Depending  on 
the  desired  transmission  rate  and  image  quality,  it  may  be  possible,  under 
the  conditions  of  moderate  motion,  to  utilize  various  psychovisual  properties 
described  in  the  previous  chapter,  such  as  spatio-temporal  response,  so  that 
the  adaptive  strategies  do  not  introduce  visible  degradations.  Nevertheless, 
for  violent  motion,  visible  degradation  may  be  unavoidable.  Therefore, 
consideration  should  be  given  as  to  what  extent  visible  degradations  can  be 
tolerated  by  a  human  viewer  for  the  application  in  mind. 

The  ideal  adaptive  strategies  should  depend  on  the  actual 
characteristics  of  the  image  signal  which  produce  the  nonuniform  information 
rate,  such  as  percentage  of  frame  motion,  velocity  of  motion,  amount  of 


spatial  detail,  etc.  However,  most  multinode  coders  use  a  single  decision 
parameter,  i.e.,  the  buffer  memory  occupancy  rate,  which  Vise  also  been  used 
under  hybrid  opcical/digital  comp ressic-n  schemes  (see  Ref.  [2]).  This  has 
lead  to  some  degree  of  success  jr.  real-time  multimode  interframe  codec 
simulations  for  coding  National  Television  System  Committee  (NTSC) 
monochrome  video  signals  at  1.5  Mbits/sec  under  minicomputer  control  [26— 28 J. 
These  can  be  called  a  class  of  multinode  coders  where  M  modes  of  coding 
operations,  such  as  were  described  in  the  previous  chapter,  are  selectively 
applied  to  a  video  signal  to  assure  its  continuous  operation,  to  provide  the 
full  available  resolution  in  the  lower  mode,  and  to  introduce  the  graceful 
and  progressive  degradation  of  image  quality  in  the  higher  mode.  The  basic 
design  problem  of  a  multimode  coder  is  lack  of  a  "genuine"  adaptive 
mechanism  to  switch  between  the  higher  and  lower  mode  of  operations. 
Instead  of  using  a  buffer  memory  occupancy  criterion,  we  propose  the  use  of 
image  activity  measures.  One  measure  of  image  activity  that  can  be  readily 
computed  by  an  optical  processor  is  image  derivative  energy  for  spatial  and 
temporal  gradients,  as  proposed  by  McCaughey  [3], 
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III. 2.  Parallel  Hybrid  Optical /Digital  Processing  Approach  To  Interframe 
Coding 

The  use  of  optical  processing  for  image  data  compression  has  been 
oriented  toward  overcoming  the  excessive  digital  processing  requirements  of 
digital  image  bandwidth  compression  techniques.  The  advantages  of  optical 
processing  of  images  are  evident  in  light  of  current  and  emerging  device 
technologies  such  as  focal  plane  CCD  or  C1D  detectors,  spatial  light 
modulators  (PROM,  LCLV',  high-speed  silicon  lithium  niobate  (SiLiNb) 
devices,  etc.),  multi-mode  optical  fibers  coupled  with  laser-based  optical 
communication  technology,  and  VLSI  technology  enabling  the  potential 
realization  of  parallel  digital  image  processors.  Within  this  new 

technological  environment,  it  is  fruitful  to  reevaluate  the  advantages  of 
both  optical  and  digital  processing  for  various  real-time  applications,  such 
as  image  data  compression  for  TV  conferencing,  military  and/or  medical 
problems  based  on  multiple  image  sequence,  etc.  The  overall  goal  is  system 
architecture  design,  i.e.,  combining  advanced  optical  and  digital  technologies 
to  improve  the  performance  of  compression  systems. 

The  viable  guideline  for  hybrid  optical/digital  system  design  would  be 
a  full  use  of  parallelism  in  image  processing.  The  conventional  and  current 
image  transmission  systems  have  been  restricted  by  a  sequential,  or  raster 
scanning  node  of  operation.  There  is  good  reason  for  it  in  terms  of  simple, 
real-time  implementation.  However,  it  may  be  useful  to  reconsider  the 
available  nodes  of  operation,  from  sequential  line-by-line  to  parallel 
block-by-block  in  spatial  domain. 

Returning  to  motion-compensated  coding  algorithms,  which  are  the  most 
sophisticated  that  have  been  developed  thus  far,  it  has  been  noted  that  they 
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suffer  L rom  excessive  search  computation  and  stability  problems.  However, 
nut  ell  the  image  subblocks  need  motion  compensation.  Also,  the  spatial  and 
temporal  gradient  information  needed  for  the  notion  estimation  can  be  easily 
cr touted  as  a  low-er.d  operation  by  the  smart' sensor  arrays  of  an  optical 
processor.  The  image  spatial  /temporal  derivative  energies  which  are  used  as 
a  criteria  for  selecting  an  appropriate  code  are  also  readily  computed  by  an 
optical  processor.  Finally,  the  hypothetical  multimode  coder,  consisting  of 
K  subimace  buffers  and  M  microprocessors,  can  be  viewed  as  a  single¬ 
instruction  multiple-data  stream/ multiple-instruction  multiple-data  stream 
(SIMD/MIKD)  parallel  proccessor.  Viewed  in  this  way,  the  adaptive  decision 
mechanism  can  be  considered  a  special-purpose  computer  (software  and/or 
hardware)  interfaced  with  both  optical  and  digital  processors,  which  takes 
information  from  the  low-end  optical  subsystem  monitoring  image  activities, 
processes  them  based  on  the  preprogrammed  decision  rules,  sends  out  the 
instructions  to  the  digital  subsystem  for  multimode  coding  as  to  which  image 
subblock  image  buffer  needs  a  particular  mode  of  processing  by  a  particular 
microprocessor,  and  generates  the  overhead  information  to  be  transmitted  to 
the  receiver. 

The  details  of  the  proposed  adaptive  multimode  predictive  coding  based 


on  parallel  image  processing  concepts  are  discussed  next. 


Hi.,';.  System  Configuration  of  Adaptive  Multimode  Coding 

Thu  basic  structure  that  we  propose  for  the  hybrid  optical/digital 
adaptive  mult inode  coder  is  given  in  Fig.  1.  The  image  sequence  is  captured 
by  a  "smart  sensor"  optical  preprocessor  that  generates  the  spatial  and 
ter.norai  gradients  for  each  subimage  block.  The  spatial  and  temporal  image 
derivative  energies  to  each  subblock  can  be  computed  by  either  a  digital  or 
an  optical  preprocessor.  Now  the  spatial  and  temporal  image  derivative 
energy  values  for  each  subblock  are  used  for  adaptive  selection  of  the 
coding  mode  applied  to  that  sunblock  or  the  temporal  derivative  and  spatial 
gradient  values  at  each  pel  may  be  used  to  implement  Pci -recursive  motion 
compensation  coding  on  a  sequential  basis  as  well. 

The  adaptive  decision  mechanism  thresholds  the  image  activity  measure 
values  for  each  block  and  generates  the  "instructions"  as  to  what  mode  of 
operation  instruction  (software  program)  is  loaded  into  which  microcrocessor 
and  which,  subimage  memory  content  is  fetched  for  execution.  When  the  Pel- 
bloc.1.:  motion  compensation  algorithm  is  applied  to  the  particular  subimage, 
the  overhead  information  about  whether  the  subimage  block  is  the  stationary 
or  moving  part  of  the  image  must  be  sent  bach  to  the  adaptive  decision 
n- onanism.  If  the  Pel-recursive  algorithm  is  applied,  the  address 
information  of  the  moving  pel  should  be  generated  on  a  sequential  basis. 
The  adaptive  decision  mechanism  can  access  the  frame  buffer  for  smoothing 
the  data  rate  in  order  to  double-check  its  performance.  Therefore,  when 
cor.  lined  with  the  overhead  information  generator,  it  is  a  solely  digital 
sub-ystem  as  is,  for  example,  a  conventional  minicomputer. 

The  parallel  image  processor  consists  basically  of  N  subblocks, 
:r..,ge/f tame  buffers/ memories  and  M  microprocessors.  Thus,  it  is  also  a 
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digital  subsystem.  M  microprocessors  for  image  processing  can  perform  P 
different  inodes  of  image  manipulations.  There  are  various  operations  that 
can  be  used  for  these  P  codes,  such  as  8-bit  PCM,  spatial/temporal 
subsaiupiing,  frame  repeating/ dropping,  motion  compensation  coding,  adaptive 
quantization,  linear  and  nonlinear  temporal  filtering,  etc.  Given  the  image 
processing  modes,  intensive  simulation  experiments  with  subjective  viewing 
should  be  used  to  find  the  appropriate  threshold  values  for  the  image 
activity  measures  that  can  be  used  for  mode  switching.  In  our  subsequent 
simulation  experiments,  five  modes  of  operation  are  used:  (1)  8-bit  PCM,  (2) 
spatial  subsanpling  (subsampling  every  two  pixels),  (3)  frame  repeating,  (4) 
motion-compensation,  and  (5)  temporal  filtering. 

Finally,  the  postprocessor  or  reconstruction  processor  at  the  receiving 
end  is  naturally  conceived  as  a  digital  subsystem  that  combines  the 
transmitted  pieces  of  subimage  blocks  and  treats  the  possible  block 
discontinuity  in  the  reconstructed  image  with  the  aid  of  the  transmitted 
overhead  information  for  subsequent  display. 
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i.  Algorithm  of  the  Proposed  Adaptive  Multimode  Coder 

Step  ].  Divide  the  image.  fields  l(X,t)  and  li(X,t+T)  into  snail 

subblocks  (say,  8X8).  Thus,  there  are  1024  8X8  subbiocks  for 
a  256*256  inage. 

Step  2.  Compute  the  temporal  differences  of  the  corresponding  ith 
subblock  between  fields  U(X,t )  and  U(X,t-fT),  such  as  L'^.C1')  » 
jlkU.t+t)  -  Ui(_X,t)j  for  i  =  1  ,...,N  =  1024. 

Seen  3.  Compute  the  temporal  image  derivative  energy  for  the  ith 
temporal  sufcblcck  by  using  the  following  formula: 


FEi  = 


d  LV.(?) 
dT 


fcr  the  i.th  subimage. 

top  4.  Compute  the  spatial  image  derivative  energy  for  the  ith 
spatial  subblock  by  using  the  following  formula: 


SE< 


S2Uj(X) 

„  n 
O^Y 


2 

for  the  ith  subimage  where  X  =  (X^,X2) 


=  iuiix1-i.>:2)  -  2ui(x1,x2;  +  uilxL+i,x2 )]2 
+  [Uif.x1,x2-i)  -  2ui(x1,x2;  +  Ui(xi,x2+i)]2. 
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Step  5. 


Step  6. 


Steo  7. 


Adaptive  decision-making  based  on  thresholding  takes  place  for 
each  subblock,  i-l,...,N  (see  Table  1).  This  rule  may  be 
changed,  depending  on  the  subjective  evaluation  tests. 

Generate  the  overhead  information  for  each  subblock  i  =  1,...,N; 
in  other  words,  information  about  which  subblock  image  i  is 
processed  by  a  particular  operation  mode  (1)  frame  repeat, 
(2)  spatial  subsampling,  (3)  motion,  compensation,  (4)  temporal 
filtering,  (5)  8-bit  PCM  ,  and  the  address  of  the  moving  and 
stationary  areas  in  the  ith  subblock,  etc. 

Reconstruct  the  block-by-block  transmitted  mosaic  image  by 
using  the  overhead  information  and,  if  necessary,  cosmetically 
treat  the  block  discontinuity  by  spatial  filtering  or  contrast 
adjustment. 


IV. 


PARALLEL  IMAGE  PROCESSING  ARCHITECTURE  FOR  MULTIMODE  CODING:  FURTHER 


RESEARCH  AREAS 

There  are  some  remaining  questions  in  the  parallel  image  processing 
approach  to  multimode  coding.  They  are  concerned  with  how  to  efficiently 
interconnect  the  subimage  buff ers/memories  with  the  multimicroprocessors. 
The  interconnection  network  problem  for  parallel  and  distributed  processing 
has  been  intensively  researched  by  the  computer  architecture  community 
during  the  past  several  years.  Siegel  et  al.  [29,30]  have  proposed  PASM:  a 
partitionable  SIKD/MIMD  system  for  image  processing  and  pattern  recognition. 
Figure  2  shows  the  block  diagram  overview  of  PASM.  There  is  a  striking 
resemblance  between  Figs.  1  and  2,  except  for  the  existence  of  the  overhead 
information  generator,  frame  buffer  and  transmission  channel.  In  other 
words,  the  parallel  digital  image  processor  in  Fig.  1  seems  to  encompass  the 
memory  storage  system,  the  memory  management  system,  the  parallel 
computation  unit,  and  the  microcontrollers  in  Fig.  2.  These  units  are 
controlled  by  the  system  control  unit,  and  correspond  to  the  digital 
decision  processor  in  Fig.  1.  Therefore,  the  same  kind  of  considerations 
toward  PASM  architecture  design  seem  to  apply  to  the  multinode  coder  design 
as  well. 

It  seems  that  the  conceptual  framework  of  parallel  digital  image 
processing  is  not  restricted  to  image  transmission,  but  can  also  apply  more 
generally  to  an  intelligent  computer  vision  application.  However,  optical 
computation  is  an  essential  part  of  preprocessing  used  in  the  computation  of 
the  spatio-temporal  gradients  and  energies  in  the  proposed  interframe 
compression  scheme.  It  has  also  been  raised  to  the  sophisticated  level  of 
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digital  computation,  provided  that  parallel  hybrid  optical/digital  systems 
architecture  is  practically  implementable. 

In  the  near  future,  the  performance  of  the  proposed  adaptive  multimode 
coding  algorithm  and  its  channel  error  effects  will  be  evaluated  by 
simulation  experiments  done  on  a  PD?  11/70  and  on  an  I~S  image  processor. 
Also,  if  time  permits,  conceptual  architecture  of  a  parallel  digital  image 
processor  for  a  multimode  coder  will  be  compared  with  PASK  architecture. 
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The  feasibility  of  optical  processing  for  ir.terframe  image  compression  has  beer-  demc.-i- 
s::i:si  oy  hybrid  optical/dicics!  ZZ?CM/"rane-or-Trana  E"C!  inter frame  compression  archi¬ 
tecture  and  its  s  inula  tier.  [1]  .  jjewever,  the  sucze=*-ve  image  quality  of  the  reconstructed 
14  th  frame  for  1.2  oits/pixel  or  1,75  bits/pixel  are  r.oo  exoeZZar.t  er.auen  to  meet  a 
broadcasting  standard.  Also,  the  transmission  rate  performances  per  second  with  the 
frame  rate:  1J  fran.es/sec  ar.d  Zoo  x  255  Zones  of  the  images  are  2.3  Mbits/ sec  or  2.4 
Mbits/sec,  respective!;',  which  are  hot  quote  competitive  acaor.se  the  1.5  Miots/sec  mor.c- 
c.-.rome  videotaZephor.a  system  peri crmar.ce  IZ].  Zr.  this  paper,  ir.  crier  to  improve  the 
hybrid  system  performance,  another  hybrid  cctical/digital  adaptive  ZTPCi/ oor.ii tocnaZ  re- 
p-er.os.hm.er.t  o.-.tsrf rams  oompressocr.  architecture  is  proposed  and  somuZated. 

3.  result  of  somuZatir.q  14  sequences  of  "alter  Crcr.hite  images  with  Zoo  x  Zae 
Zones  ar.o  Z  oit  intensify,  reconstructed  images  w o th  “r.  excellent  image  quality  are 
obtained  with  1.2*  -its/sec  fer  the  frame  rate:  20  frames/sec,  ar.d  the  average  ccnpres- 
socr.  ratic  per  formsr.ee  is  15. 34.  u,v*'t  ri^  p.m  v-jiti  i*  c  7  b  c>i  . 

Zr.troducticn 

Zr.age  data  cemprassicr.  is  a  very  active  topic  cf  research  ir.  image  processing  [2’., 

■;4'i .  Traditionally,  mest  image  data  ccmpressocr.  schemes  are  dcmo.-.atad  by  digital 
processing  due  tc  the  fact  that  digital  systems  are  inherently  flexible  ar.d  reliable, 
"ewever,  it  has  beer,  demcr-s orated  that  optical  processing  car.  oe  used  fer  image  data 
compression ,  oy  Hunt,  et.al.  Col,  [oj,  and  [*j.  One  cf  the  schemas  for  image  data 
compression  architecture  t.-ich  cptical  im.plemer.taticr.  is  the  ZT?CZ:  method  [5].  this  is  a 
data  compression  scheme  which  functions  analogously  tc  conventional,  digital  OPCZI  compres¬ 
sion,  except  t.-.at  t.-.e  specific  compression  steps  are  implemented  by  incoherent  optical 
processing.  The  extension  of  the  ZZPCZS  method  into  ir.terframe  compression  using  a  frar.e- 
to-frame  digital  ZPC:  scheme  leads  ti  hybrid  opticai/digital  architecture  with  electrc- 
cptica-  devices  ,  [31,  such  as  PF.OM  or  CCS  or  liquid  crystal  light  value. 

-he  bcot.ler.ecx  of  this  proposed  hybrid  cpoicai/digitai  ir.terframe  compression  scheme 
is  mainly  caused  oy  the  motion  displacement  contained  ir.  neighboring  frames  as  well  as 
the  frame  prediction  mechanism,  the  predicted  frame  used  in  [7]  is  ; use  a  previous  frame, 
assuming  the  motion  displacement  bec.ver.  fram.es  is  small.  This  is  .not  usually  me  case 
ir.  real  ir.carfrare  image  sequences.  Therefore,  motion  displacement  errors  between  neigh¬ 
boring  frames  are  accimulated  as  time  ar.d  frar.es  go  or.,  ar.d  unacceptable  distortion  is” 
generated  m  the  reconstructed  image  if  the  lew  frequency  quantitation  level  is  not  large 
enough,  to  allow  the  low  frequency  component  wnicn  constitutes  most  cf  me  metier,  displace¬ 
ment  to  ce  correctly  coded.  Also,  era  surpective  ir.3ge  quality  cf  the  reconstructed  14th 
frame  for  1.2  cits,  pixel  cr  1.75  bits/pixe.  is  net  excellent  or. tug.-,  to  meet  a  broadcasting 
star.iard.  Zr.  addition,  t.-.e  transmission  rate  performances  ter  second  with  the  frame  rate: 

ZZ  frames,  sec  ar.d  256  pixels  x  256  lines  cf  the  images  are  Z .  3  .Icits/sec  cr  2.4  .Ibits/sec, 
respectively.  This  is  net  quite  competitive  against  the  1.5  .‘ifcits,  sec  of  monochrome  video¬ 
telephone  system  performance'  [2;.  Z.-.  this  paper,  in  order  to  improve  the  hybrid  optical/ 

digital  system  performance ,  another  hybrid  optical/digita.  adaptive  ZTPCM/ conditional 
replenis.br.ent  in oa r frame  oompressicr.  architecture  is  proposed  ar.d  simulated.  T.-.e  conven¬ 
tional  conditional  rapler.is.tnsr.t  oacnn.que  is  based  or.  t.-.e  segmentation  of  the  moving 
pixels  from  t.-.e  stationary  pixels,  and  .only  t.-.e  significantly  moving  pixels  are  transmitted 
along  with,  their  address  information .  .At  t.-.e  racai— ..-.g  and,  t.-.e  stored  prsoeedir.g  image 
is  repeated  with  t.-.e  significantly  moving  pise.s  replaced  oy  t.-.e  transmitted  moving  tixa.s. 
Zha  following  are  problems  frequently  .sscciated  with  t.-.e  conditional  replenishment  tac.v- 
•.true:  1,  ouffer  overflow,  ar.d  Z  address  coding.  They  will  not  oe’ discussed  m 

datai.,  except  t.-.at  the  t.-.res..r.i  Ze-  e  .  cf  t.-.e  metier,  ieteotcr  is  recursive set  sc  t.-.at 
the  buffer  car.  store  t.-.e  .artsst  number  cf  moving  pixels  wit.-.m  me  buffer  canacitv. 


T.-.e  essentia-  tart  if  t.-.e  cor.r 
there,  car.  be  fividad  into  direct 
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sprite*”  car.  be  incorporated  ir.  the  conditional  replenishment  technique  by  using  various 
modes.  The  above  motion  detection  ar.i  adaptive  strategies  will  re  discussed  ir  the  del  - 
Icwmg  secticr  under  the  ccrstrairt  cf  the  hybrid  optical, 'digital  implementation . 

Ir  the  Simla  tier,  experiment  ,  we  examine  the  performances  of  four  systems  including 
the  proposed  adaptive  hybrid  cptical/di vital  II?  Ol/cor.diticnal  rep  ler.  i  s.-ment  system. 

£  vs  ter.  ccr.f  icuratior.  of  a  temtcrallv  adattive  hybrid 
cotiral  cicital  ir.terfreme  cc.T.rressicr.  scher.e 

ZTi -S reel 2  - 1  ve  D?C!  sririil  co~cr?ss i c o  sci^srr.s 

Ir.terfrar.e  coding  is  the  rest  effirier.t  image  coding  syster.  ir  terr.s  of  the  irate  data 
reduction  rate  because  spatial  redundancy  as  well  as  tsr.pcral  redundancy  car  be  eliminated 
if  the  reverent  involved  between  neicnbcri.ro  frames  is  assumed  to  be  small.  As  or.e  of  the 
effioiert  spatial  compress ior.  systems,  it  has  been  demonstrated  that  the  interpolate ve 
I?C‘.  ;  lire./  method  can  rarcrstruco  an  image  wich  excellent  quality  "i:il  «=  0.137%)  with 
a  bit  rate  cf  '  bits/pixel,  which  is  better  than  a  conventional  IPO)  method,  figure  1(a) 
shews  the  incoherent  optical  system  implementation  of  IIPCI)  compression  and  figure  Kb) 
shows  the  digital  implementation  of  lirCf.  compression.  The  main  differences  between  IlrC, 
and  l?i::  are  as  follows. 

(1)  ll?c:;  does  not  assume  causality  of  a  pixel  sequence  in  an  image  line,  that  is, 
ire:  is  a  pixel-to-pixel  operation  as  opposed  to  II?C)  which  takes  a  parallel 
image  plane  to-- image  plane  operation. 

(1)  the  pixel  estimate  of  I?d  is  replaced  with  a  neighborhood  smoothing  of  the 
entire  image  p.a.-.e,  which  is  called  a  low  frequency  image  in  IIPC.. 

(3)  the  pixel  difference  cf  TPC)  is  replaced  with  the  difference  between  the 
original  image  and  the  estimated  ;r  smoothed  low  frequency  image,  which 
is  called  a  hig.t  frequency  image  in  TTPC1. 

1-.)  another  significant  difference  is  a  non-  feedback  structure  in  ITPCt  as  opposed, 
to  a  feedback  structure  in  TPC).  The  feedback  cf  a  pixel  estimate  in  .TPC*.  is 
used  to  cancel  the  accumulated  quantitation  error  from  the  quantized  difference. 
However,  by  using  two  channel  incoherent  optical  syster.  implementation,  the 
feedback  structure  car.  be  eliminated  without  accumulating  the  quantization 
error.  This  is  cue  to  the  fact  that  the  visibility  cf  the  distortion  at  the 
receiver  is  reduced  by  the  two-channel  structure  in  which  the  high  frequency  image 
containing  various  distortion  is  added  to  the  low  frequency  image. 

’-he  similarities  between  "PC:  and  TPC)  can  be  summarized  in  the  following  two  points 
which  lead  to  a  common  data  compression  strategy 

(1!  the  low  frequencies  m  the  image  car.  be  predicted,  removed  from,  transmission 
by  a  difference  step,  and  then  repeated  at  t.-.e  receiver. 

’I)  the  high  frequencies,  which  are  r.cb  eliminated  by  prediction  ar.d  differencing, 
pccsess  only  a  small  amount  of  the  total  image  energy  ar.d  car.  be  accurately 
coded  with  a  small  rummer  cf  bits. 

Ir.  Hybrid  Optical.  Tigital  Implementation,  the  sue  sampling  ar.d  interpolation  cperat.cn 
cir.  -a  raalicai  ar.d  zr.e  lif isrer.cir.c  zzer aricr.  car.  ir.-lerr.er.tad  dv  ar.  2 -■rrc-r  i 2te 

The  spatial  compression  performances  or  the  IIPC:  system,  for  every  Ir.d  pixel  ar.d  every 
•' t-h  pixel  3ur sampling  ar.d  interpolation  are  giver,  in  table  1(a)  ar.d  lib),  when  1TVQ  ar.d 
71  levels  are  changed  from  1  to  e'4.  The  overall  optimal  performance  (M'.SI  ar.d  bit  rate  per 
ixel)  seems  tc  re  giver,  by  the  every  Ir.d  pixel  subsamplir.g  and  interpolation  with  ITl'l  » 

.'..-I  •  -  ant  .1)11  «  1.1371:%.  This  combination  cf  quantization  level  parameters  has 
been  used  for  tr.s  subsequent  adaptive  ir.tsrdram.e  compression  simulation. 


5 3 1 c r.  scr.emes 


:mprassicr.  s o.-.eme  ar.d  simulated  ar.d  the  results  are  giver,  in  .me  resu.es 

:cw  that  unless  t.-.e  ITT  level  is  large  enough  IT 1  »  3  or  1  = ;  ,  the  accumulated  rrt.cr 
i.splacemer.  t  errtrs  among  a  sequence  cf  frames  causes  severe  distortion  ar.d  art.  farts 
t.-.e  reconstruct'd  image  rf  the  1-th  frame.  Also  t.-.e  sui;active  image  rual.t-  :: 
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broadcasting  standards,  and  the  transmission  rate  per  second  w 
and  256  pixels  :<  156  lines  cf  sequence  ad  the  inaces  is  1.3*3. 
fair  image  quality. 


.h  33  frames/sec  d r are  rate 
hbits/sec  with  reasonably 


It  has  been  observed  that  the  videcccr.derence  images  involve  mainly  snail  movement  so 
that  the  stationary  pixels  cf  the  preceeding  frame  car.  be  repeated  to  reconstruct  the  next 
drar.e  at  the  receiver,  and  only  the  moving  pixels  are  trar.sr.itt <*d.  in  other  words,  the 
.conditional  reoler.ishr.ent  codir.c  can  be  exoected.ta  perdorr.  better  than  Frame -to-Frane 
CPC1  coding  in  terr.s  cd  both  KI-iSE  ar.c  transmission  bate. 


The  typical  syster.  components  cd  the  conditional  replenishment  coding  system  ere 
(1)  notion  detector  or  3igr.idicant  change  detector,  (2!  pixel  selector  ,  based  on  the 
address  information  produced  by  (1)  ,  (3)  unidorr,  quantizer  (64  levels)  ,  (4)  address 

coder,  (5)  drar.e  budder,  end  (6)  reconstructer. 


The  essential  part  cd  the  conditional  replenishment  coder  is  the  motion  detector  or 
.significant  change  detector,  which  can  be,  in  general,  divided  into  two  categories:  direct 
and  indirect  methods.  The  direct  method  is  based  on  laborious  calculation  of  the  velocity 
[3!  ,[13]  or  the  estimation  o d  motion  displacement  cd  a  moving  pixel  between  neighboring 
frames,  using  a  3  dimensional  image  or  dield  model  [11],  [12-]. These  methods  are  the  bases 
for  notion  estinatior./cor.pensation  coding  whose  signidieantly  improved  perdcrmar.ce  has  been 
demonstrated  dor  spatially  ur.idcm  translated  displacement  by  Xetravali,  et.el  [13]  and 
ethers.  2r.  order  to  implement  conditional  replenishment  interdrame  ceding  combined  with 
:TPC:  spatial  compression  within  the  framework  od  the  hybrid  oetieal/digital  system,  an 
.ndirect  method  od  detecting  signidieant  change  between  neighboring  dramas  by  threshold¬ 
ing  the  drame  didderer.ee  i3  predered  to  the  direct  method  ir.  which  the  velocity  and/cr 
motion  displacement  are  calculated  from  the  pixel-to-pix«I  operation.  Also,  the  state  of 
-  budder  to  smooth  moving  pixel  data  generated  by  movement  car.  be  used  as  an  indirect 
r.dicatisr.  od  the  significant  change.  .In  addition,  NMSE  od  the  frame  difference  can  in- 
.icate  movement  activity. 


"ithin  tha  framework  of  hybrid  optical/digital  architec 
e  designed  as  a  hybrid  system,  in  which  two  liquid  crystal 
iffarence  image  ar.d  the  fecal  plane  charge  injection  device 
image  intensity  into  a  digitally  addressable  signal.  Then, 
the  adaptive  scheme  car.  be  made  by  using  a  microprocessor  ac 


cure,  the  motion  detector  can 
light  valves  generate  a  drame 
(CIO)  converts  the  difference 
the  decision  making  needed  for 
ting  on  the  digital  signal 


res  the  CIO. 


The  adaptive  strategies  can  be  incorporated  immediately  into  conditional  replenishment 
coding,  which  is,  indeed,  a  partially  adaptive  scheme  itself,  by  accommodating  several 
- idderer.t  modes.  They  are  as  follows: 

(1)  Adaptive  change  od  the  threshold  level  depending  on  the  movement  activity 
and  the  buddsr  state. 

(2)  Seduction  od  spatial  and/or  temporal  resolution  by  using  subsamplir.g  or 
filtering. 

(3)  Suspension  cd  replenishment  due  to  large  movement. 

(4)  Field  repeating  due  to  very  small  or  no  movement. 

3iver.  the  constraint  od  the  optically  implemented  spatial  ccmpressior.  system,  the 
^beve  adaptive  strategy  (2)  has  a  limited  utility  because  UPC:  is  not  quite  suitable  dor 
-.ar.ging  the  spatial  resolution  od  the  limited  region  (moving  area)  in  the  image.  Also, 

.  :e  image  plane- to- image  plane  operation  od  I3FC1  makes  it  inefficient  to  directly  calculate 
the  velocity  or  moticr.  displacement  od  signidieantly  moving  pixels  as  in  digital  processing. 
Therefore,  we  adopted  an  indirect  change  detection  method  using  adaptive  thresholding 
!  the  drame  difference  based  on  the  budder  state  and  the  adaptive  strategy  (3)  the 
isper.sien  od  the  replenishment,  to  dorr,  adaptive  hybrid  opticai/digital  I3PC!/ccniitior.al 
replenishment  ir.terdram.e  compression  architecturs  as  shown  in  Figure  2. 


Simulation  experiments  a.-.d  results 


To  dsr.onstrats  the  superior  performance  od  the  adaptive  hybrid 
conditional  replenishment  interdrame  compression  architecture,  a  se 
•r?«sl"sa-*  were  carried  cut  in  the  digital  Image  Analysis  lacerate 
#  «rsr*4  • 


sttical/digitai  IDPCl/ 
les  cd  digital  sinsiiatien 
;  cd  the  University  od 


The  source  data  dor  the  simulation  consisted  cf  a  sequence 
from  a  television  broadcast  od  Waiter  Crcr.kite.  The  frames  are 
i  xeis  resolution,  with  e  bits  of  intensity  per  pixel. 


cf  14  digitised  frar.ee 
digitised  et  236  x  236 
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Assuming  that  address  coding  is  accurate/  no  channel  coding  has  been  used,  and  no  error 
has  been  added  to  the  transmission  channel,  we  have  simulated  the  following  four  systems 
for  eer.pariscr.  of  the  compression  performances  (bits/see,  with  frame  rate  of  30  frar.es/sec) 
and  image  quality 

(A)  A  r.onadaotive  conditional  replenishment  system  with  the  1st  frame  transmitted 
by  64  level  PC!. 

(B)  A  nonadaptive  hybrid  optical/digital  ISPCI/cenditienai  reolenishment  system 
with  the  1st  frame  transmitted  by  lire:  (KTS  -  1TTQ  -  32). 

(C)  An  adaptive  conditional  replenishment  system  with  64  level  PCI. 

(C)  An  adaptive  hybrid  optieal/digital  ISPCt/ccnditional  replenishment  svstem 
(lira  -  2,  LTCQ  -  32)  . 

The  adaptive  strategies  are  based  or.  frame  replenishment  suspension  or  frame  repeating 
depending  or.  the  :3iS2  of  the  successive  frame  difference.  In  other  words,  if  XISE  of  the 
frame  difference  is  less  than  0.05%,  the  frame  is  repeated,  and  if  IZiSZ  is  more  than  0.5%, 
the  frame  rsolensihment  is  suspended  and  onlv  a  spatial!'/  compressed  image  by  64  level  PCI 
or  ICPCt  (iri'Q  -  32,  K~Q  -  2)  is  transmitted*. 

the  motion  detection  strategy  consists  of  frame  differencing  ar.d  iterative  threshold¬ 
ing  as  shewn  in  rigure  3,  whose  initial  threshold  is  set  as  the  mean  value  of  the  frame 
difference  and  cr.e  hundredth  of  the  variance  of  the  frame  difference  is  incremated  till  the 
maximum  numbed  of  moving  pixels  are  detected  within  the  maximum  buffer  sice  (3030  pixels) 
for  uniform  quantisation. 

Table  2  summarised  the  transmission  rate  (bits/sec) ,  the  average  bit  rate  per  frame 
ibits/frar.e)  and  the  average  compression  per  frame  ratio.,  of  four  systems.  Table  3  shows 
:~:CZ.  the  number  of  moving  pixels,  transmitted,  compression  ratio,  and  the  cumber  of  bits 
transmitted  for  the  proposed  adaptive  hybrid  optieal/digital  IBPCi/conditional  replenish¬ 
ment  system. 

Sraph  1  shows  how  the  NhSl  performances  of  four  systems  change  along  14  frames. 


as  in  the  ease  of  Hybrid  optieal/digital  I3PC:/rrams”to--rrane  3PCI  ceding  jystem,  the 
motion  displacement-  between  neighboring  frames  seems  to  be  accumulated.  This  could  be  due 
to  the  fact  that  the  buffer  size  is  limited  so  that  sufficient  number  of  moving  pixels  ere 
not  correctly  detected  and  transmitted  to  yield  the  reliable  reconstruction  at  the 
receiving  end. 

2r.  the  adaptive  scheme,  in  the  same  way  as  above,  I3PC ’  based  conditional  replenish- 


cor.strueted  pictures  by  both  systems  are  excellent  as  the  14th  frame  reconstruction  is 
shown  in, Picture  2.1,  and  2.2.  The  worst  reconstructed  images  ere  the  6th  frames  shown  i- 
Fi-.urs  i.3^  ant  -.cr.parir.g  t.-.e  original  6t.-.  frams  and  14th  frame  shown  in  Picture  J.l. 

3...  with  the  reconstructed  ones,  it  is  observed  that  a  so-called  ghost  effect  slightly 
legraded  the  facial  part  but  not  the  background  of  the  reconstructed  image. 

It  is  r. otewerthy  that  the  frame  difference  sometimes  contains  noise  in  the  stationary 
area  as  shown  in  Picture  4 . 1  This  r.oise  can  be  filtered  out  in  order  to  reduce 

:he  number  of  moving  pixels  to  be  buffered. 

Conclusions  ar.d  further  research 

have  describee'  ar.  ir.terfr&me  data  centre ssior.  svstem  which  uses  an  adaptive 
treeegy  of  replenishment  suspension  and  2  bi'ts/pixel  ICPC:  spatial  compression  scheme 


itnin  the  framework  of  Hybrid  optieal/digital  architecture.  The  results  of  this  temporal 
-y  adaptive  mterframe  ccmpressicr.  simulation  are  encouraging  because  of  the  overall  good 
ob;ective  and  subjective  image  quality  CZlZZ  -  3.08  *  1. £T%) ,  the  average  number  of  bits 
•r  f-'-  Kbits/frame ..  the  average  compression  ratio  per  frame.  16.84  ,  and  the  trans- 

issior.  rata  with  frame  rate  33  fram.es/sac;  1.322  :  lb  its . 


There  are  seme  obvious  directions  in  which  future  research  in  this  scheme  should  b 
i.-.ted.  "or  example,  the  motion  detection  scheme  can  be  modified  by  filtering  the  lew 
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contrast  moving  areas  ar.d  the  frame-to-frame  noise  in  the  static  areas  o f  the  picture,  or 
by  an  adaptive  reduction  in  the  sensitivity  of  the  notion  detector.  The  next  step  of  rais¬ 
ing  the  efficiency  and  adaptability  of  the  interfrar.e  cor.pressier.  system  seems  to  be  the 
development  of  the  more  efficient  and  real-time  motion  detection  processor  within  the  frame 
work  of  hybrid  cptical/digital  architecture. 

Also,  it  is  './orth  examining  the  minimum  sice  of  the  frame  buffer  in  order  that  the 
subjective  image  quality  of  the  image  reconstructed  by  the  r.onacaptive  conditional  replen¬ 
ishment  coding  is  sufficiently  good.  Usually,  the  buffer  sice  is  constrained  by  the  pro¬ 
gressing  digital  memory  technology  and  cost  so  that  it  may  cease  to  be  the  major  limiting 
factor  to  design  the  interfrar.e  coder. 

Tinally,  it  is  important  to  incorporate  the  effect  of  address  coding  and  channel  coding 
under  the  channel  error. 
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Picture  2.2 
adaptive  IDPC-i  based 
conditional  replenishner.t 
14th  frame  reconstruction 


Picture  2.3  Picture  2.4 

adaptive  PCI  besed  adaptive  13 PCI  based 

conditional  replenishner.t  conditional  replenishment 
6th  frane  reconstruction  6th  drace  reconstruction 
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Picture  3.2 
:rlpir.sl  14th  drama 


Picture  3.1 
original  Ss 


Picture  4 . 1 
frame  difference 
between  6th  and  7th  drar.es 
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Radiometric  and  Spatial  Transformations 
for  Image  Data  Compression 

Research  Objectives 

Adaptive  image  data  compression  has  been  exclusively  a  digital 
processor  function  because  of  the  nonlinear  and  space- vari ant  opera¬ 
tions.  However,  a  scanning  optical  system  can  achieve  geometric 
space- variant  operations  by  alteration  of  the  optical  system’s 
instantaneous  field  of  view  (IFOV)  and/or  variation  of  the  path  of 
the  IFOV.  Radiometric  s pace- vari ance  can  be  achieved  by  the  usage 
of  a  programmable  mask,  e.g.,  a  PROM  or  liquid  crystal.  In  the 
research  under  this  task  we  are  exploring  the  extent  to  which  radio- 
metric  and  geometric  space- variance  can  improve  the  performance  of 
a  simple  non-adaptive  optical  compression  scheme,  IDPCM,  previously 
developed  under  sponsorship  of  Grant  AF0SR- 76- 3024 . 

Appendix  Contents: 

This  appendix  consists  of  a  summary  of  research  results  on 


employing  radiometric  and  spatial  transforms,  plus  a  paper  on  the 
theoretical  basis  for  deriving  spatial  transformations. 
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1)  INTRODUCTION 

An  optical  means  for  intraframe  compression,  IDPCM,  originally 
developed  by  Hunt1,  is  modified  by  pre-compression  and  post-compression 
processing.  The  pre-compres s i on  processing  takes  the  form  of  a  rad¬ 
iometric  and/or  spatial  transformation  performed  on  the  original  image 
to  produce  an  image  with  stationary  mean  and  energy  (radiometric)  and/ 
or  autocorrelation  (spatial),  respectively,  as  defined  over  sub-blocks 
of  the  image.  The  post-compression  processing  is  simply  the  inverse  of 
whatever  transformation  that  was  applied  to  the  original  image. 

Such  a  modification  to  the  basic  compression  step  is  possibly 
a  way  to  lower  bit  rates  for  given  image  quality,  or  improve  image 
quality  at  a  given  bit  rate. 

The  analysis  that  follows  describes  the  effect  of  various  rad¬ 
iometric  transformations  on  objective  and  subjective  image  quality, 
for  different  compression  bit  rates.  A  spatial  transformation  is  also 
investigated. 
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2)  PROCEDURE 

The  image  used  in  this  analysis  is  the  128  X  128  pixel,  256  grey 
level  (8  bit/pixel)  Walter  picture  shown  in  Fig.  1.  The  small  size  is 
necessary  to  keep  processing  time  to  a  minimum.  A  256  X  256  8  bit/pixel 
image  is  also  processed  for  a  few  of  the  cases  described  below,  to  de¬ 
termine  if  image  size  is  important. 

Radi ometri c 

The  general  procedure  is  as  follows;  the  original  image  is  com¬ 
pressed  with  the  IDPCM  architecture  alone  for  various  values  of  low 
frequency  quantization  value  L  and  high  frequency  quantization  value 
H,  to  define  a  set  of  images  that  are  purely  the  result  of  the  bit  com¬ 
pression.  Then  another  set  of  images  are  derived  for  the  same  L,H  val¬ 
ues  as  in  the  first  set,  but  with  the  further  processing  steps  of 
a)  performing  a  radiometric  transformation  on  the  image  before  the  com¬ 
pression  step  and  b)  performing  the  inverse  radiometric  transforma¬ 
tion  on  the  compressed  image  from  part  a.  The  radiometric  transform  used 

.  2 

is  that  derived  by  Strickland  ,  namely: 

g( j  ,k)  =  -s [f ( J ,k)  -  yN]  +  vs 
n 

di 


where  f(j,k)  is  the  original  value  of  the  jth,kth  pixel, 

gU.k)  is  the  transformed  value  of  the  jth,  kth  pixel, 
v  N  is  the  original  value  of  the  mean  determined  over  each 

sub-block  of  the  image  lin  all  transformations  described 
below,  the  sub-blocks  are  8X8  pixel  squares;  there  are 
256  of  these  sub-blocks  in  the  128  X  128  image  used! 


y  is  the  new  value  of  the  sub-block  mean  (the  same  for 
al 1  sub-bl ocks ) 

oN  is  the  original  standard  deviation  determined  over  each 
sub-block,  and 

is  the  new  value  of  the  sub-block  standard  deviation 
(the  same  for  all  sub-blocks). 

The  subscripts  s  and  N  refer  to  stationary  and  non-stationary , 
respectively.  By  stationary,  it  is  meant  that  the  statistics  (namely 
mean  and  standard  deviation)  of  each  sub-block  after  the  above  trans¬ 
formation  will  be  the  same  as  those  of  any  other  sub-block.  The  par¬ 
ameters  of  this  transformation,  ys  and  o  ,  are  these  stationary  values; 
they  are  set  by  the  user.  (The  parameter  os  is  actually  determined  by 
setting  the  desired  stationary  energy,  Rs;  i.e., 


Rs  -  yj 


(2) 


-  thus  it  is  the  energy  and  mean  parameters  that  are  discussed  in 
what  f ol 1 ows . ) 

This  transformation  is  applied  to  each  sub-block  individually; 
a  blocky  image  would  result  if  no  further  modification  was  made. 
Therefore  the  coefficients  and  (one  set  for  each  block)  are 
linearly  interpolated  between  block  centers  to  yield  a  smooth  image. 
An  example  of  this  radiometric  transformation  acting  on  Fig.  1  is 
shown  in  Fig.  3b,  with  y  =  100  and  R  =  12000.  The  histograms  of  the 

V  J 

original  and  transformed  images  are  shown  in  Figs.  6  and  8,  respect- 
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Table  1  lists  the  various  reconstructed  images.  For  given  values 
of  L  and  H  the  energy,  R$,  of  the  radiometric  transformation  is  var¬ 
ied  to  determine  what  effect  the  transformation  h-.s  on  the  compression 
process.  The  transform/compression/inverse  transform  process  for  Fig. 

3a  is  illustrated  by  the  sequence  of  figures  outlined  in  Table  2.  The 
inverse  transformation  is  simply  accomplished  by  inverting  eq.  (1); 
i.e.,  solve  for  f(j,k).  Values  of  L  from  2  to  16  are  looked  at  to  see 
the  sensitivity  of  the  low  frequency  quantization  to  the  transforma¬ 
tion.  Note  that  in  almost  all  of  the  transformations,  the  stationary 
mean  value  is  100.  This  value  is  chosen  arbitrarily;  the  primary  in¬ 
terest  here  is  to  study  the  effect  of  a  given  mean,  varying  the  energy 
parameter . 

The  high  frequency  quantization  value,  H,  was  chosen  arbitrar¬ 
ily  at  first  to  match  the  low  frequency  quantization  value,  L,  as  shown 
in  Table  1.  A  subsequent  analysis  (see  Figs.  19  and  20)  is  performed 
to  isolate  the  effect  of  the  high  frequency  value  H  from  the  compres¬ 
sion  and  low  frequency  quantization  processes,  and  is  discussed  in  the 
results  section. 

Spatial 

The  discussion  above  refers  to  the  radiometric  transformation. 

2 

A  spatial  transformation,  derived  by  Strickland  ,  was  applied  to  the 
compression  process  in  the  same  type  of  way  as  the  radiometric  trans¬ 
form;  i.e.,  spatially  transform  the  original  image/compress  the  result/ 
inverse  transform  the  compressed  image.  The  spatial  transformation  has 
the  property  of  normalizing  the  autocorrelation  width  of  each  sub¬ 
block  ever  the  entire  image,  where  again  the  sub-block  is  an  8  X  8 
pixel  square.  The  shape  of  the  autocorrelation  function  for  each  block 
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is  approximated  by  an  elliptically  symmetric  function.  The  width  of 
this  function  is  estimated  for  each  block  in  the  vertical  and  horiz¬ 
ontal  directions;  from  this  information  each  sub-block  of  the  image 
is  in  effect  expanded  to  varying  degrees  in  the  horizontal  and  verti¬ 
cal  directions,  depending  on  the  original  autocorrelation  of  the  sub¬ 
block  in  the  two  orthogonal  directions.  Thus  sub-blocks  with  high  cor¬ 
relation  (uniform  grey  areas)  are  left  alone,  and  sub-blocks  with  low 
correlation  (edges)  are  expanded.  Fig.  16  is  an  example  of  the  spatial 
transform  used  in  this  analysis.  Fig.  17  uses  the  spatial  transform 
alone  with  the  compression  step,  and  Fig.  18  uses  both  the  spatial  and 
a  radiometric  transformation  with  the  compression  step. 

NMSE 

A  quantitative  comparison  of  two  images  is  the  normalized  mean 
square  error,  or  NMSE.  It  is  defined  as  follows: 

t  (fr(j,k)  -  f(f,k))* 

NMSE  =  ^ - 5 - 

i  f(j.kr  (3) 

j  ,k 


where  frlj,k)  is  the  value  of  the  jth,  kth  pixel  of  the  reconstructed 
image,  and  f(j,k)  is  that  value  for  the  original  image.  The  NMSE 
values  are  given  in  Table  1  for  each  reconstructed  image,  where 
Fig.  1  was  used  as  the  original  image. 
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3)  RESULTS;  RADIOMETRIC 

L  =  2,  H  =  2;  Low  frequency  channel;  128  X  128  Walter  image 

The  L  =  2,  H  *  2  case  will  be  investigated  first,  with  regard  to 
the  radiometric  transformation.  Looking  at  the  compressed  image  with 
L  *  2 ,  H  =  2  (Fig.  2)  with  no  radiometric  transform  applied,  the  low 
frequency  quantization  is  readily  apparent  by  the  two  level  pixel  val¬ 
ues;  bright  pixels  correspond  to  pixel  values  in  the  original  image 
that  were  above  some  threshold,  dark  pixels  to  those  that  were  below 
the  threshold.  Looking  at  Fig.  3a  (radiometric  transform  with  energy 
=  12000  used),  a  similar  quantization  of  pixel  values  is  obvious,  again 
due  to  the  threshold  imposed.  (The  intermediate  frames.  Figs.  3b  to 
3g,  used  in  the  production  of  Fig.  3a  will  be  referred  to  in  the  sub- 
suquent  analysis-  see  Table  2.)  However,  in  the  latter  image,  the  IDPCM 
induced  bright  and  dark  areas  (quantization  noise)  are  uniformly  spread 
over  the  image.  This  is  the  effect  of  the  stationary  transformation  on 
the  original  picture  block  means,  and  the  two-level  quantization.  The 
amount  of  this  low  frequency  quantization  noise  can  be  seen,  looking 
at  the  progression  of  Figs.  4, 3a, and  5,  to  be  controllable  by  adjust¬ 
ing  the  parameter  of  energy,  R$,  from  10,500  to  12.U00,  to  20,000,  res¬ 
pectively.  This  parameter  is  ultimately  affecting  the  histogram  width 
of  the  transformed  image  (before  compression),  and  this  width  af¬ 
fects  the  position  of  the  uniform,  equally-spaced  pixel  bins  in  the 
uniform  quantization  step  (low  frequency  channel).  This  may  in  fact 
be  the  most  important  effect  of  the  energy  parameter  for  low  quantiza¬ 
tion  values.  For  the  case  here,  two  quantization  levels,  a  slight  change 
ir.  the  upDer  or  lower  limits  of  the  histogram  to  be  quantized  causes 

tue  breakpoint  pixel  value  to  change.  It  so  happens  that  in  this  image. 


the  pixel  values  of  the  background  wall  are  very  close  to  this  cutoff 
value,  and  their  quantized  representation  is  thus  sensitive  to  par¬ 
ticular  histogram  upper  and  lower  bounds,  which,  as  stated  above,  are 
related  to  the  energy  parameter  of  the  transform.  (The  greater  the 
energy,  the  greater  the  spread  in  upper  and  lower  bounds.)  Note  that 
Fig.  4,  with  its  different  statistics,  has  a  breakpoint  value  below 
the  wall  pixel  values,  and  the  wall  is  uniformly  bright.  See  Fig.  7 
and  Fig.  8  for  examples  of  stationary  image  histograms  (before  com¬ 
pression)  of  Figs.  4  and  3a,  respectively  (the  original  image  histo¬ 
gram  is  in  Fig.  6). 

Another  effect  of  the  radiometric  transform  is  visible.  Recall 
that  Fig.  3g  is  the  reconstructed  low  frequency  frame  before  the  inverse 
transform.  Applying  the  inverse  transform  to  this  image  added  to  3f  produces 
Fig.  3a,  the  final  frame;  note  the  grey  levels  introduced  into  Fig.  3a. 

Dark  areas  are  made  brighter,  and  vice  versa.  (This  is  the  inverse 
stationary  process  acting  on  the  block  means.) 

The  effects  of  the  radiometric  transformation  for  the  case  1=2, 

H  =  2  in  the  low  frequency  channel  are  now  summarized: 

1)  Dark  areas  of  the  picture  normally  lost  in  the  compression  step 
are  brightened  (and  vice  versa),  thus  the  low  frequency  quant¬ 
ization  noise  is  redistributed  more  uniformly  over  the  image. 

2)  The  value  of  the  quantization  breakpoint  depends  upon  upper 
and  lower  bounds  of  the  histogram,  which  is  related  to  the  mean 
and  energy  parameters  of  the  transformation;  notice  the  dif¬ 
ference  between  Figs.  4  and  3a. 

3)  The  inverse  transform  re-introduces  a  range  of  grey  levels  that 
are  lost  in  the  quantization  process. 


The  above  discussion  refers  to  the  low  frequency  channel;  the 
high  frequency  channel  is  now  discussed,  with  respect  to  the  single 
frame.  Fig.  3a. 

Fig.  3f  is  the  Laplace  quantized  Fig.  3e;  these  are  intermedi¬ 
ate  frames  that  correspond  to  Fig.  3a.  Fig.  3f  is  the  image  added  to 
Fig.  3g,  the  interpolated  low  frequency  image,  to  give  the  image  which 
is  then  inverse  transformed  to  produce  Fig.  3a.  Note  that  the  quant¬ 
ized  image  is  a  poor  representation  of  the  high  frequency  image:  the 
significant  edges  are  present  in  the  quantized  image,  but  the  added 
quantization  noise  is  uniform  over  the  picture,  degrading  the  informa¬ 
tion  content.  Thus,  in  all  the  frames  shown  at  H  =  2  (L  doesn't  af¬ 
fect  high  frequency  quantization),  the  high  frequency  channel  contri¬ 
butes  little  information  to  the  final  picture.  Note  that  in  the  L  =  2, 

H  *  2  sequence  land  in  other  sequences),  the  high  frequency  noise  is 
visible  in  the  background. 

To  determine  the  number  of  quantization  levels  needed  to  render 
an  accurate  high  frequency  image,  H  was  varied  from  2  to  4  to  6  to  8, 
in  the  compression  of  an  untransformed  picture  that  was  very  similar  to 
Fig.  3e.  The  significant  results  are  in  Figs.  19  (H=4),  and  20  ( H*  8 ) . 
The  subjective  fidelity  of  the  quantized  high  frequency  image  becomes 
acceptable  at  H  =  b  (this  image,  not  shown,  is  almost  identical  to 
Fig.  20).  The  effect  of  low  high  frequency  quantization  levels  should 
be  kept  in  mind  in  the  discussion  of  the  subsequent  frames  (when  H  is 
less  than  6). 

Quantitative  comparison 


The  normalized  mean  square  error,  or  NMSE,  is  the  quantitative 


no 


comparison  of  the  original  image  with  the  reconstructed  image.  The 
NMSE  of  Fig.  2  and  Fig.  4  are  .131  and  .069,  respectively;  an  improve¬ 
ment  of  50"  is  apparent  (see  Table  1).  As  the  energy  parameter  in¬ 
creases  above  10,500  (see  Figs.  3a,  5),  so  does  the  NMSE.  The  reason 
the  NMSE  is  lower  in  the  transformed  images  is  evident  by  noting  that, 
in  the  untransformed  image,  much  of  the  wall  area  of  the  picture  has 
been  discarded  due  to  the  low  frequency  quantization  step;  however, 
the  transformed  image  has  a  more  unifc;mly  represented  wall  (depending 
on  the  particular  pixel  distribution  that  determines  the  quantization 
breakpoint).  Also,  the  inverse  transform  adds  grey  levels  to  the  image 
(corresponding  to  the  grey  levels  in  the  original),  which  tends  to 
improve  the  NMSE. 

L  =  4,  H  =  4 

Again  increasing  the  energy  from  10,500  (Fig.  9)  to  12,000  (Fig. 

10)  shows  the  alteration  of  the  image  that  is  characteristic  of  the 
shift  in  the  breakpoint  values.  Now,  however,  four  levels  of  quant¬ 
ization,  thus,  three  breakpoint  values,  are  being  used  in  the  quant¬ 
ization.  The  value  of  the  center  breakpoint  is  close  to  the  single 
breakpoint  derived  for  the  two  level  quantization,  thus  the  wall  pixel 
pattern  resembles  that  of  Fig.  3a.  Qualitatively,  the  quantized  nature 
of  the  untransformed  image.  Fig.  11,  has  been  removed  in  Fig.  9.  Quant¬ 
itatively,  there  is  a  corresponding  reduction  in  the  NMSE  of  the  two 
images,  from  .024  to  .016,  respectively. 

L  ■  8,  H  *  8 

Here  again,  the  quantized  nature  of  the  untransformed  Image,  Fig. 
12,  is  smoothed  by  the  radiometric  transform  producing  Fig.  13;  however, 


the  NMSE  shows  virtually  no  change,  form  .004  to  .003,  respectively. 


L  -  lfe,  H  ■  16 

There  is  virtually  no  change  in  using  the  radiometric  transfor¬ 
mation  for  such  a  high  low  frequency  quantization  (L=16),  as  seen  by 
comparing  Fig.  14  ( untransformed)  and  Fig.  15  (transformed).  The  NMSE 
is  the  same  (.005)  for  both  cases.  The  low  value  for  H  (=2)  indicates 
that  all  of  the  information  present  is  from  the  low  frequency  chan¬ 
nel;  the  high  frequency  channel  noise  can  be  seen  superposed  on  the 
image  for  both  cases.  It  is  interesting  to  see  how  much  information 
is  present  in  the  low  frequency  channel.  Subsampling  at  a  coarser  rate 
may  be  worthwile  in  reducing  the  bit  rate  for  such  high  values  of  L. 

256  X  256  Walter  image 

The  128  X  128  Walter  image  was  used  in  the  above  sequences,  to 
reduce  computing  time.  Two  256  X  256  images  were  studied  with  L  =  2, 

H  =  16,  to  check  for  consistency,  and  gave  results  similar  to  those  of 
the  128  X  128  images  described  above. 

Spatial  transform 

Preliminary  results  of  the  spatial  transform  shown  in  Fig.  16  for 
the  l  *  4,  H  *  4  sequence  are  Figs.  17  and  18.  Fig.  17  is  the  result  of 
the  spatial  transform  alone;  Fig.  18  is  a  combination  of  the  spatial 
and  a  radiometric  transform.  The  results  show  little  effect  due  to  the 
spatial  transform  used;  a  more  radical  spatial  transform  is  probably 
requi red . 
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4)  CONCLUSIONS 

A.  Radiometric  transformations  seem  to  work  best  (qualitatively) 
to  improve  compressed  frame  quality  for  low  frequency  quantiza¬ 
tions  of  4  and  8  (i.e.,  L*4,  L*8). 

B.  Breakpoint  values  in  the  uniform  low  frequency  quantization  step 
are  sensitive  to  image  statistics  (namely,  lowest  and  highest 
pixel  values),  and  thus  change  with  different  radiometric  trans¬ 
formations.  Such  changes  cause  dramatic  (image  dependent)  changes 
in  the  final  reconstruction;  e.g.,  the  wall  pixels  in  the  Walter 
image . 

C.  Quantization  noise  introduced  by  the  compression  is  more  uni¬ 
formly  distributed  over  the  image  when  the  radiometric  transform 
is  applied. 

D.  The  inverse  radiometric  transform  re-introduces  a  range  of  grey 
levels  that  are  lost  in  the  quantization  process. 

E.  The  high  frequency  Laplace  quantized  Image  contributes  little  to 
the  final  reconstructed  picture  for  H  less  than  6. 

F.  The  spatial  transform  (warp)  tested  did  little  to  change  the 


qualitative  value  of  the  compressed  image. 
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5)  FUTURE  WORK 

A.  Improve  low  frequency  quantization  algorithm;  uniform  quant¬ 
ization  appears  too  crude.  The  radiometric  transform  may  produce 
a  histogram  whose  low  frequency  quantization  can  be  optimized 
(e.g.,  Laplacian,  Gaussian).  If  not,  can  the  radiometric  trans¬ 
form  be  taylored  to  produce  an  optimum  histogram  shape? 

B.  Try  subsampling  at  a  coarser  rate  in  the  low  frequency  channel 
when  the  low  frequency  quantization  number  is  large;  say  L  greater 
than  or  equal  to  16. 

C.  Try  more  radical  spatial  transformations  to  provide  a  better  test, 
for  this  type  of  transform;  also  determine  the  effects  of  coars¬ 
er  sampling  rates  used  in  conjunction  with  spatial  transforms. 

D.  Add  spatially  uniform  random  noise  to  the  warped  image  before 
the  compression  step,  subtracting  the  exact  same  noise  after  the 
compression;  this  will  reduce  contouring,  at  the  price  of  a  gran¬ 
ular  image.  Applying  the  inverse  warp  may  then  reduce  the  effect 
of  the  granularity  in  image  areas  effected  most  by  the  warp  (i.e., 
areas  of  low  correlation). 


E.  Combine  spatial  and  radiometric  transformations. 


TABLE  1  :  List  of  Figures 
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Figure 

kill 

Mean 

Energy 

Descri pti on 

NMSE 

1 

— 

— 

original 

— 

2 

2,2 

— 

IDPCM  alone 

.131 

3a 

2,2 

100 

12000 

IDPCM, radiometric 

.090 

4 

2,2 

100 

10500 

IDPCM, radiometric 

.069 

5 

2,2 

100 

20000 

IDPCM, radiometric 

.119 

6 

— 

— 

histogram  of  Fig.  1 

— 

7 

— 

— 

histogram  of  Fig.  4 

— 

8 

— 

— 

histogram  of  Fig.  3a 

— 

9 

4,4 

100 

10500 

IDPCM, radiometric 

.016 

10 

4,4 

100 

12000 

IDPCM,  radi ometri c 

.023 

11 

4,4 

— 

IDPCM  alone 

.024 

12 

8,8 

— 

IDPCM  alone 

.004 

13 

8,8 

100 

11000 

I DPCM ,radi ometri c 

.003 

14 

16,2 

— 

IDPCM  alone 

.005 

15 

16,2 

127 

43000 

IDPCM, radiometric 

.005 

16 

— 

— 

spatial  only 

-  -  -  - 

17 

4,4 

— 

IDPCM, spatial 

.039 

18 

4,4 

100 

10500 

IDPCM, rad. .spatial 

.019 

19 

-.4  . 

— 

high  freq.  quant. 

— 

20 

-.8 

— 

high  freq,  quant. 

TABLE  2  :  Intermediate  Processing  Steps  of  Figure  3a 
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Step 

1 

2 

3 

4 

5 

6 

7 

8 
9 


Figure 

1 

3b 

3c 

3d 

3e 


3f 

3g 

not  illustrated 
3a 


Description 
origi nal 

radiometrically  transformed  Fig.  1 

(mean  =  100,  energy  =  12000) 

subsampled  Fig.  3b 

uniformly  quantized  Fig.  3c 

high  frequency  of  Fig.  3b  (derived  from 

subtracting  the  linearly  interpolated  Fig. 

3c  -not  illustrated-  from  Fig.  3b) 

Laplace  quantization  of  Fig.  3e 
linearly  interpolated  Fig.  3d 
the  sum  of  Fig.  3f  and  Fig.  3g 
the  inverse  transform  of  the  above  image; 
the  final  image 
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1.  INTRODUCTION 


The  statistical  behavior  of  images  is  of  fun¬ 
damental  importance  m  many  areas  of  image 
processing.  Incorporation  of  a  priori  statis¬ 
tical  knowledge  of  the  spatial  correlation  of 
an  image,  for  example,  can  lead  to  consider¬ 
able  improvement  in  many  image  restoration 
algorithms.  The  recent  trend  in  nonstation- 
ary  restoration  techniques  (1),;2)  has  been 
paralleled  by  a  similar  shift  in  data  compres¬ 
sion  schemes,  from  spatially  nonadaptive 
differential  pulse  code  modulation  (DPCM)  sys¬ 
tems  to  adaptive  DPCM  systems  where  the  coef¬ 
ficients  of  the  optimum  predictor  are  tuned  to 
the  local  statistical  behavior  of  the  image. 
The  superiority  of  spatially  adaptive  schemes 
is  now  well  supported  (3). 

The  foregoing  evidence  in  support  of  using 
nonstationary  statistical  image  models  leads 
us  to  pose  the  following:  suppose  we  turn 
the  problem  around— spatially  transform  or 
warp  an  image  so  that  it  possesses  stationary 
statistics  before  feeding  it  to  a  nonadaptive 
process ,  such  as  nonadaptive  DPCM  or  the  non¬ 
adaptive  Wiener  filter — and  then  use  the 
reverse  warp  transformation  to  restore  the 
original  image  geometry.  We  might  then  ex¬ 
pect  to  see  improved  performance  from  those 
image  processes  that  assume  stationarity . 

In  a  previous  paper.  Hunt  (4)  proposed  the 
use  of  nonstationary  statistical  image  models 
in  data  compression.  Schemes  were  outlined 
for  the  measurement  of  space-varying  image 
parameters  and  implementation  of  image  warp¬ 
ing  using  hybrid  digital/optical  hardware. 

The  present  paper  represent!  a  continued 
effort  towerd  developing  trsnsforme  for  gen¬ 
erating  statistically  stationary  imagas.  Our 
emphasis  will  bs  on  the  applicstion  of  such 
transforms  for  enhancing  the  performance  of 
DPCM  date  compression,  although  we  also  in¬ 
tend  to  investigate  their  gains  in  imaga 
restoration  applications.  Throughout  tbs 
paper  we  use  digital  simulation  to  tast  the 
transformations,  although  we  anticipate  thee 
ultimately  they  can  be  implemented  optically. 

2.  CONDITIONS  FOR  STATIONARY  STATISTICS 

We  are  concerned  with  generating  image*  pos¬ 
sessing  wide-tense  stationarity,  involving 
only  first-  and  second-order  moment*  (S). 

For  a  two-dimensional  image  f(x,y),  these  ere 
the  mean  and  autocorrelation  statistic!, 
respectively 


for  stationarity  of  the  process  f(x,y)  in  the 
wide-sens*  are  aa  follows: 
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»  f(j+r,k*s)  »  Rjjlr,s).(4) 


In  other  words,  the  means  of  the  neighborhoods 
N  are  constant,  and  thsir  autocorrelation 
functions  have  constant  lag  parameters . _  M  is 
the  dimension  of  a  square  neighborhood.'  In 
Eq.  (4) ,  an  autocorrelation  matrix  with  dimen¬ 
sions  (2M-l,Jt)  is  sufficient  due  to  redun¬ 
dancy  inherent  in  the  reel  and  even  function. 
Thi*  property  further  allows  us  to  specify 
the  autocorrelation  by  four  unique  profiles 
shown  in  Figurs  1  and  designated  aa  NS,  NE, 

Ew,  and  SE  as  points  on  the  compass. 

The  spatial  variance  of  the  image  statistics 
is  reflected  by  the  spatial  variance  of  the 
mean  and  autocorrelation.  In  the  ca«*  of  the 
mean,  we  have  a  single  number  for  each  neigh¬ 
borhood.  Somewhat  more  complex  is  the  auto¬ 
correlation  which,  being  e  function  of  space, 
is  charaetsrized  by  its  energy,  width,  and 
shapa.  Henca,  a  total  of  four  parameters  may 
be  used  to  describe  neighborhood  statistics. 
According  to  our  specifications  for  station¬ 
arity  given  in  Eqs.  (3)  and  (4) ,  any  varia¬ 
tion  between  the  perimeters  for  all  neighbor¬ 
hoods  will  result  in  nonstationary  behavior. 

3.  THEORY  OF  TRANSFORMATIONS  FOR 

- BTWTgfnav~gfA«5TtCS - 

3.1.  Maan  and  Msan-Squars  Energy 

Consider  the  ease  of  a  squara  neighborhood  N 
of  f ! j ,k)  consisting  of  M*  pixels.  We  assume 
that  .y  and  Ry (0, 0)  are  initially  nonststion- 
ary  with  respect  to  other  neighborhoods.  In 
general,  we  anticipate  that  stationary  sta¬ 
tistics  will  be  generated  by  a  simple  linear 
transformation  in  intensity  of  the  form 

g  ( j  ,k)  -  Aflj.k)  ♦  B.  (5) 


-j.(x.y)  •  f(x,y)  dxdy 

N 

Rylx.y,;,-)  *  ^  f (x,yl f (x*C ,/♦-) 

N 

■  dxdy , 
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Her*  the  subscript  N  denotas  that  tha  moment 
statistics  are  measured  in  local  nsighbor- 
hoods.  For  tha  digital  cast,  the  conditions 


Ths  conditions  for  stationarity  in  g(j,k)  are 

expressed  as 
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where  _s  and  R«(0,2!  are  stationary  values. 
Solving  Eqs .  (5)  and  (6)  for  constants  A  and 
S  yields. 


v  -  v°-o>r  .  f. 

.-N;  *  V°'0,J  'S 


(7a) 


stationary  standard  deviation 
*  nonstationary  standard  deviation 


B  *  -j  “  Ais.  (7b) 

Equations  (5)  and  (7)  combine  to  give  the 
transformation  for  stationary  mean  and  mean- 
square  energy,  or  stationary  mean  and  vari¬ 
ance: 

?<:,*)  *  r  Cf(j.Jc)  -  a  ]  *  w_  (8) 

N 


Equation  (8)  represents  a  filtering  algorithm 
that  is  identical  to  an  algorithm  developed 
for  the  spatial  contrast  enhancement  of 
images  ■'  6 )  . 


It.  practice,  Eq.  >8)  will  generate  images 
containing  intensities  that  fall  outside  the 
available  dynamic  range  of  the  original  image, 
in  our  case  B  bits,  or  gray  levels  0  to  2S5. 
Negative  data  are  obviously  impractical  for 
optical  implementation.  Furthermore,  we  wish 
to  restrict  the  upper  gray  level  limit  to  255 
so  as  not  to  contradict  our  ultimate  aim  of 
data  compression.  We  will  therefore  redefine 
the  transformation  in  Eq.  (8)  to  accommodate 
the  dynamic  range  constraint:  0  s  g(J,k>  s 
255.  Tne  new  transformation  is 


?<:.*>  •  *[f<j.it>  -  uN]  -  .s, 

where 

*  *  r  if  0  :  g(j,k)  i  255 


otherwi se 

k  •  minimum 


255 
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(9a) 

(9b! 


( 9c) 


where  [f !  J  ,k> -JN]max  *"<»  C * <  j  ,k> -JN3mln  are 
the  maximum  positive-going  and  maximum 
negative-going  variations  (about  the  mean) , 
respectively.  By  definition.  Eq.  (9!  will 
generate  images  with  stationary  mean  and 
approximately  stationary  mean-square  energy. 


3.2.  Autocorrelation  Shape 

In  section  X  we  characterised  the  shape  of 
the  autocorrelation  by  four  unique  profiles. 
Clearly,  transformation  to  stationary  behav¬ 
ior  requires  that  we  assign  some  kind  of  model 
to  each  profile  or.  alternatively,  a  model 
for  the  autocorrelation  ae  a  whole.  An  image 
field  is  considered  to  be  a  good  fit  to  a 
first-order  Markov  process  for  which  the  cor¬ 
relation  between  individual  pixels  is  propor¬ 
tional  to  their  geometric  separation  (3).  we 
choose  to  neglect  the  diagonal  profiles  in 
Figure  1  and  define  an  elllptically  symmetric 
form  of  the  autocorrelation. 


R_w  «  F.J..C,-)  »  constant  '  expf^j^,  t  ,  ] 

(10) 

where  R»js  and  Rrw  are  the  profiles  in  orthog¬ 
onal  North-South  (vertical)  and  East-West 
(horizontal)  directions  respectively. 

Stationary  autocorrelation  width  along  a  par¬ 
ticular  direction  is  produced  by  normalizing 
each  nonstationary  : ,  measured  in  tne  same 
direction,  to  some  value  c*.  It  is  logical 
to  equate  with  the  minimum  value  of  z, 
corresponding  to  the  highest  correlation. 

The  correlation  widths  are  subsequently  nor¬ 
malized  by  a  resampling  (interpolation)  of  the 
image  data  by  a  factor  £,  given  by 


s 


E«j  is  therefore  the  expansion  factor  (along  a 
given  direction  in  the  image)  required  to 
transform  the  autocorrelation  width  t>(j  of  a 
neighborhood  N  into  a  stationary  value  cg.  The 
two-parameter  model  given  by  Eqs.  (10)  gives 
rise  to  expansion  along  orthogonal  directions: 
Ejjs  and  Egw-  We  derive  methods  for  perform¬ 
ing  local  spatial  expansions  in  section  5. 


4 .  APPLYING  TRANSFORMATIONS 
FOR  atATIOKARY  S’fATiSf ttS 


4.1.  Test  Image 


Figure  2  shows  the  256  «  256  pixel  8-bit 
"Walter”  image  on  which  we  base  our  results. 
The  statistics  are  visibly  nonstationary. 

For  example,  examining  contiguous  blocks  of 
16  *  16  pixels,  we  find  that  the  means  range 
from  13  to  203,  and  the  mean-squares  range 
from  172  to  43386.  The  global  mean  is  90, 
with  dynamic  range  from  9  to  249.  For  most 
of  the  results  to  follow  we  will  consider 
neighborhoods  of  16  «  16  pixels,  giving  a  total 
of  256  contiguous  blocks  throughout  the  test 
scene.  The  resolution  of  local  statistical 
behavior  is  not  optimum  at  this  block  size 
since,  since  ideally  we  would  like  to  resolve 
single  edges  in  a  block. 

4.2,  Mean  and  Mean-Square  Energy 

we  recall  that  the  basic  transformation  to 
produce  stationary  mean  and  mean-square  en¬ 
ergy  is  given  by  Eq.  (8).  Applying  dynamic 
range  constraints  gave  a  practical  transform, 
Eq.  (9).  Before  applying  this  we  need  to  say 
something  about  the  effects  of  noise.  Let  us 
assume  for  algebraic  simplicity  that  any  noise 
present  is  additive,  with  zero  mean.  We  can 
write  the  noisy  image  as 

f' (j,k>  -  f (j,k)  *  n(j,k) .  (12) 

which,  when  substituted  for  f(j,k)  in  Eq.  (8) 
gives 

9<J.*>  -  ^  Cf(j.k)  *  n(j,k)  -  u„]  *  u#. 

(13) 

From  Eq.  (13)  we  can  see  that  the  noise  term 
is  amplified  by  the  factor  3a/sp.  Ideally, 
we  want  the  noise  to  remain  in  the  background. 
In  practice,  we  can  partially  attain  this  goal 
by  rewriting  Eq.  (13)  as  follows: 

g(j,k)  •  K[r  (j,k>  -  aN]  *  a,  (14) 


RjjS  •  RjjC.Oi  •  constant  •  exp[-:NS:  ;  ]  where 
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uses  so  generate  two  values  of  expansion  for 
each  neighaorhhod,  as  follows: 


Constant  T  is  a  threshold  related  to  the  rms 
noise  level  _-R.  Assuming  tnat  f'j.k)  and 
are  uncorreiated  random  variables, 
then  setting  T  »  .-n.I  would  be  equivalent  tc 
saying:  if  the  rms  signal  variations  are 
less  or  equal  to  the  estimated  noise,  do  not 
attempt  to  equalize  the  mean-square  energy, 
or  variance.  Instead,  adjust  only  the  mean. 
The  same  threshold  condition  on  K  can  also 
oe  applied  to  Eos.  i 9 ) .  giving: 
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(Note:  the  neighborhood  subscript  has  beer, 

omitted  from  Eq.  (16)  for  the  sake  of  clar¬ 
ity.)  Figures  4a  and  b  are  maps  of  E^s  and 
Eew  derived  from  the  transformed  image  in 
Figure  2,  in  which  each  value  of  expansion  is 
represented  by  a  block  of  16  «  16  pixels  with 
uniform  gray  level.  The  data  are  scaled  to 
3  bits ,  corresponding  to  the  range  of  expan¬ 
sion  values. 


5.  IMAGE  WARPING  EASED  ON  EXPANSION  MAPS 


we  now  address  the  problem  of  performing  geo¬ 
metrical  warping  to  accommodate  each  one  of 
256  local  expansions.  A  general  approach  is 
to  assume  an  m-order  polynomial  model  for  the 
spatial  distortion,  giving. 
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when  r  jj  -  T  and 
K  *  1 

when  -  v  i  T. 

Figure  3  shows  the  test  image  after  trans¬ 
formation  using  Eq.  (15)  with  fit  •  127  and 
: s  *  165.  equivalent  to  a  stationary  mean- 
square  of  433B6  (as  original  image) .  We  use 
:r  »  5  as  cur  noise  estimate.  The  final 
image  would  have  exhibited  discontinuities 
at  the  borders  between  neighborhoods  if  we 
had  used  single  values  of  and  cN  for  each 
block  of  16  «  16  pixels;  instead,  we  use  bi¬ 
linear  interpolation  between  neighborhoods 
to  provide  unique  values  for  every  pixel. 

The  result  is  a  continuous  image.  Alterna¬ 
tively,  we  could  have  measured  .(j  and  jjj  in 
neighborhoods  centered  on  all  256  »  256  pixels 
ir.  the  test  image,  obviously  at  the  expense 
of  increased  computation.  In  Figure  3,  the 
mean  is  constant  at  127  (measured  in  16  «  16- 
pixel  blocks) ,  and  the  mean-square  energy 
ranges  from  16,300  to  26.000.  We  may  con¬ 
clude  that  the  image  is  stationary  in  these 
statistics  to  a  good  approximation.  The 
deviations  in  the  actual  mean- squares  about 
the  intended  value  of  43386  are  due  to  the 
noise  thresholding  and  dynamic  range  con¬ 
straints  . 


in  m-  •  , 

q  »  l  l  bMxly  <  (17b) 

i»0  j«0 

where  (p.q)  and  (x.y)  are  coordinates  of  the 
two  spaces  (warped  and  unwarped,  or  vice 
versa  since  the  transformation  can  take  place 
in  either  direction) . 

A  method  for  deriving  the  polynomials  is  pro¬ 
vided  by  superposing  the  expansion  of  each 
neighborhood  individually,  with  the  constraint 
that  the  influence  of  each  expansion  on  the 
relative  geometry  of  other  neighborhoods 
should  be  minimized.  The  technique  we  have 
developed  is  based  on  control  point  grids. 

5.1.  Control  Point  Happing 

Methods  of  using  control  points  in  similar 
geometrical  rectification  and  warping  appli¬ 
cations  are  well  documented  (7,8).  In  our 
case,  we  define  a  regular  grid  of  control 
points  whose  coordinates  coincide  with  the 
comers  of  each  square  neighborhood,  giving, 
for  the  previous  example,  a  total  of  17  »  17 
control  points.  Our  problem  is  to  determine 
how  this  regular  grid  in  the  original  image 
space  is  redistributed  by  the  complex  inter¬ 
actions  of  local  expansions  centered  on  each 
neighborhood.  As  mentioned  above,  we  begin 
by  determining  how  the  grid  is  affected  by 
a  single  expansion  applied  to  the  first 
neighborhood.  We  then  apply  a  second  expan¬ 
sion  to  the  second  neighborhood,  and  so  on. 


4.3.  Expansion  Maps  for  Stationary 
Autocorrelation  Shape 

In  section  3.2  we  discussed  the  principles  of 
equalizing  autocorrelation  shape  in  neighbor¬ 
hoods  by  performing  localized  direction- 
dependent  expansions.  After  transformation 
to  stationary  mean  and  variance,  the  auto¬ 
correlation  functions  of  contiguous  16  •  16- 
pixel  neighborhoods  are  computed,  least- 
squares  fit  values  of  and  : ew  are  com¬ 
puted,  corresponding  to  the  autocorrelstion 
profiles  in  orthogonsl  directions  (defined 
in  Figure  1!  ■  The  minimum  values  of  for 
each  direction  from  all  neighborhoods  are 


Each  single  expansion  has  two  components  EjjS 
and  Eew.  acting  in  orthogonal  directions, 
shown  in  Figure  5.  The  region  to  be  expanded 
is  contained  within  a  circle  of  radius  d/.'l. 
Expansion  causes  vector  displacements  at  a 
point  (i)  lying  on  the  circumference  of  this 
cirele  a>  cording  to  the  following: 
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(18) 

Outside  the  circle,  for  example  at  point  (ii) , 
the  transformations  in  x  and  y  are 

p  »  —  eosS-IE^-l)  *  x  (19a) 
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Wb  see  iron  Eos.  i!8)  that  the  four  neighbor- 
r.ood  control  points  lying  within  the  circle 
are  expanded  by  Ess  in  y  end  Epk  ir.  x.  The 
two-dimensional  transformation  does  not  pre- 


serve  the  original  geometry 


the  remaminc 


neighborhoods .  We  can  see  intuitively  that 
to  do  so  is  impossible.  Ir.  the  following 
section  we  describe  a  technique  for  reducing 
the  mutual  interference  of  local  neighborhood 
expansions. 

5. 2.  Iterative  Control  Point  Mapping 


We  have  extended  the  preceding  method  of  suc¬ 
cessive  expansions  to  incorporate  correction 
factors  to  offset  the  affects  of  nonlocalized 
expansions.  In  Figure  6  we  define  the  geo¬ 
metrical  terminology  for  neighborhoods  of 
arsitrarv  shape.  The  expansion  components 
in  this  scheme  are  defined  as 


The  iterative  algorithm  for  control  point 
mapping  then  proceeds  as  follows: 

Notation: 

Neighborhoods:  N.  ,N: , .  . .  ,N  ;,6 
Oesired 

expansions:  (E^-.E^l  j  ,  (E  -,£  )  2 ,. .. 

Measured  NS  ew  ns  EW 

expansions :  (E^s  .E^)  ,  ,  lE^.E^)  2  , . . . , 

(ENS'EEWf  :S6 

Algorithm: 

Step  1  ,  , 

Measure  (E„C,E_,)  ;  for  N .  . 

Step  2  NS 

Apply  corrected  expansions  E^.  /E^_  , 
eSW./eEW.  t0  »»•  S:  “NS, 


Step  511 

Measure  (Ens.Eew)256  *or  N:56 
Step  512 

Apply  expansions  ENS  ;5s/ENSi5S. 

EEW;:j/EeWI56  to  N;55 

Return  to  Step  1,  repeat  for  I  iterations. 

In  the  above  we  are  applying  correction  fac¬ 
tors  Ens , Egw  to  the  desired  expansions  E^s, 
Eew-  The  effectiveness  of  the  algorithm  is 
measured  by  the  convergence  of  the  rms  ex¬ 
pansion  errors  measured  at  the  end  of  each 
iteration.  In  Figure  T  we  have  the  warped 
control  point  map  for  our  previous  two-com¬ 
ponent  expansion  data  after  nine  iterations. 
Figure  8  shows  the  converging  rms  error 
curves.  There  is  apparently  little  to  be 
gained  by  continuing  beyond  nine  iterations. 

5.3.  Control  Grid  Interpolation 

At  this  point  we  have  two  sets  of  control 
points:  the  first,  a  regular  grid  in  origi¬ 
nal  unwarped  image  space:  the  second,  a 
severely  warped  grid  that  will  ultimately 
establish  the  geometry  of  a  stationary  image. 
The  quadrilateral  vertices  leontrol  points) 
map  directly  to  the  corresponding  square 
vertices.  Interior  points  are  mapped  accord¬ 
ing  to  interpolation  based  upon  the  mapping 
of  the  vertices.  We  use  local  mapping, 
based  on  small  sets  of  neighboring  control 
points.  For  example,  we  can  express  three- 


poir.t  mapping  between  triangles  as 
p  «  a.  *  a.x  *  a-.y , 

(21) 

q  »  br  -  b;X  *  b;y 

Ir.  Eqs.  (21)  we  create  two  sets  of  three  lin¬ 
ear  equations  in  the  same  number  of  unknowns 
by  specifying  that  the  vertices  of  a  triangle 
map  into  the  vertices  of  the  corresponding 
figure  in  warped  space.  Solution  of  these 
equations  leads  to  the  mapping  coefficients 
a  and  b.  Global  warping  based  on  two  sets 
of  17  «  17  control  points  is  therefore  achieved 
by  computing  32  »  22  polynomial  transforma¬ 
tions  that  are  applied  individually  to  corre¬ 
sponding  neighborhoods. 

Although  up  until  now  we  have  talked  in  terms 
of  (p,q)  being  warped  space  and  (x,y)  being 
unwarped  space,  in  actual  fact  we  derive  the 
transformations  after  reversing  the  defini¬ 
tions.  Hence,  for  each  integral  coordinate 
in  (x,y)  in  output  warped  space,  the  trans¬ 
formations  tell  us  where  to  look  in  unwarped 
(p,q)  space.  If,  as  is  the  general  ease, 
the  coordinates  (p,q)  are  nonintegral,  then 
we  have  to  resort  to  gray  level  interpolation 
to  generate  an  output  pixel.  Bilinear  inter¬ 
polation  of  the  four  nearest  neighbors  is 
satisfactory.  By  defining  the  direction  of 
the  transformation  from  warped  to  unwarped 
space,  we  make  the  task  of  pixel  interpolation 
easier. 

Figure  5  shows  the  warped  image  derived  by 
three-point  spatial  trnas formation  of  the 
intermediate  image  in  Figure  3,  based  on  the 
control  point  map  in  Figure  7.  It  is  inter¬ 
esting  to  note  that  noise-smoothing  by  low- 
pass  filtering  would  now  cause  relatively 
less  blurring  of  visually  important  edges  due 
to  the  spatial  expansion  of  such  regions. 

This  is  obviously  a  simple  example  of  the  ad¬ 
vantages  of  processing  stationary  images. 

6.  CONCLUSIONS 

We  have  discussed  transformations  for  produc¬ 
ing  images  with  wide-sense  stationary  first- 
»nd  second-order  moment  statistics.  These 
transforms  are  reasonably  straightforward  tw 
derive  and  to  apply.  We  note  that  the  in¬ 
verse  spatial  transformation  is  simple  to 
accomplish  since  it  involves  only  switching 
the  roles  of  the  two  sets  of  control  points 
when  implementing  Eq.  (21)  .  Optical  imple¬ 
mentation  of  forward  and  reverse  spatial 
transformations  in  intensity  in  Eq.  (15)  by 
optical  means  is  less  obvious;  it  may  be 
necessary  to  resort  to  hybrid  digital/optical 
technology.  Further  work  will  be  directed 
at  testing  the  effectiveness  of  stationary 
images  in  DP CM  data  compression,  compared  with 
nonstationary  images. 
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Figure  1  Autocorrelation  function  profile 
designations. 
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Fioure  2  Original  test  image. 


(bl  E-W  data 

Figure  4  Expansion  ratios  map:  elliptically 
symmetric  autocorrelation  model 


Figure  5  Principle*  of  two-dimen*ional 
control-point  napping. 


Figure  7  iterative  warping:  control  point 
grid  after  nine  Iteration*  (ex* 
pen* ion*  of  Figure  S) 


